---
title: "Radical Weakness --- Do Radical Parties Receive Fewer Ministries?"
subtitle: "Replication File"
author:
    -   name: Timo Sprang
        orcid: 0009-0005-3798-7747
        email: sprang@soz.uni-frankfurt.de
        affiliation:
        -   name: Goethe-Universität
            city: Frankfurt am Main
            region: DE
format:
    html: 
        filters: 
            - parse-latex
fig-dpi: 300
---

# Main Article

```{r}
#| label: setup

library(dplyr)
library(ggplot2)
library(rnaturalearth)
library(marginaleffects)
library(ggpubr)
library(forcats)
library(sampleSelection)
library(tibble)
library(modelsummary)
library(lmtest)
library(sandwich)
library(knitr)
library(kableExtra)
library(broman)
library(tidyr)
library(ggdark)

options(modelsummary_factory_html='kableExtra')
options(modelsummary_factory_markdown='kableExtra')
set.seed(2111)

if (knitr::is_html_output()){
    theme_set(dark_theme_bw())
    lc = "black"
    ltop = NULL
    hl = FALSE
    css = "border-bottom: 1px solid"
} else if (knitr::is_latex_output()){
    theme_set(theme_bw())
    lc = "black"
    ltop = "scale_down"
    hl = TRUE
    css = NULL
}
```

```{r}
#| label: functions

change_coefs <- function(pseudo, model, vars, type, cluster = df$country_name){

    cluster <- coeftest(model, vcovCL, cluster = cluster)

    if(type == "sel"){
        pseudo <- summary(pseudo)
        pseudo$coefficients[,1] <- summary(model)$estimate[vars, 1]
        pseudo$coefficients[,2] <- cluster[vars, 2]
        pseudo$coefficients[,4] <- cluster[vars, 4]
    } else if(type == "out"){
        pseudo <- summary(pseudo)
        pseudo$coefficients[,1] <- summary(model)$estimate[c(vars, nrow(summary(model)$estimate)), 1]
        pseudo$coefficients[,2] <- cluster[c(vars, nrow(summary(model)$estimate)), 2]
        pseudo$coefficients[,4] <- cluster[c(vars, nrow(summary(model)$estimate)), 4]
    }

    pseudo
}

make_matrix <- function(data, nvar, colnames){
    x  <- matrix(rnorm(nrow(data) * nvar), ncol = nvar) |>
        data.frame()
    colnames(x) <- colnames

    x
}

heck_gof <- function(model, first, last){
    gofs <- model |>
        glance() |>
        t() |>
        data.frame() |>
        rownames_to_column() |>
        rename(v1 = 2) |>
        filter(rowname %in% c("AIC", "logLik")) |>
        mutate(
            v2 = "",
            rowname = ifelse(rowname == "AIC", "AIC", "Log Likelihood")
        ) |>
        select(rowname, v2, v1) |>
        rbind(list("Country FEs", "Yes", NA))
        

    attr(gofs, "position") <- c(first:last)

    gofs
}

imr_sel <- function(model){
    (dnorm(predict(model, part = "outcome", type = "unconditional")) / pnorm(predict(model, part = "outcome", type = "unconditional")))
}
```

```{r}
#| label: data-descriptives

df <- read.csv("df.csv")
dfc <- df |>
    filter(cabinet_party == 1)

party_obs <- nrow(df)
coa_obs <- n_distinct(df$cabinet_id, na.rm = T)
country_obs <- n_distinct(df$country_name_short, na.rm = T)
```

```{r}
#| label: method-infos

maxyear <- max(df$year)
country_list <- gsub(
    ", United",
    ", and the United",
    gsub(
        "Netherlands",
        "the Netherlands", 
        paste(sort(unique(df$country_name)), collapse = ", ")
    )
)

relabs_cor <- format(round(cor.test(df$parlmean_dif, df$left_right)[[4]], 3), nsmall = 3)
relabs_p <- scales::pvalue(cor.test(df$parlmean_dif, df$left_right)[[3]], add_p = TRUE)
relabs_n <- format(nrow(df), big.mark = ",")

relabs_cor_g <- format(round(cor.test(dfc$parlmean_dif, dfc$left_right)[[4]], 3), nsmall = 3)
relabs_p_g <- scales::pvalue(cor.test(dfc$parlmean_dif, dfc$left_right)[[3]], add_p = TRUE)
relabs_n_g <- format(nrow(dfc), big.mark = ",")
```


```{r}
#| label: fig-example-mean
#| fig-height: 3
#| fig-width: 7
#| fig-cap: "Relative radicalism. **Note**: Hypothetical parties on the left-right axis, depicting the operationalisation of the relative radicalism measure."
#| fig-pos: t

data.frame(
    type = c(rep("Party", 4), "Mean"),
    party = c("A", "B", "C", "D", "Mean"),
    lr = c(2, 5, 7, 8, ((45 * 2) + (5 * 40) + (7 * 10) + (8 * 5)) / 100),
    seats = c(45, 40, 10, 5, 100),
    distance = c(2, 1, 3, 4, NA)
) |>
mutate(
    label1 = ifelse(
        type == "Party",
        paste0(party, "\nSeats: ", seats),
        "Centroid"
    ),
    label2 = ifelse(
        type == "Party",
        paste0("LR: ", lr, "\nDist.: ", distance),
        paste0("LR: ", lr)
    )
) |>
ggplot(aes(x = lr, y = 0, colour = party, shape = type, size = seats)) +
    geom_hline(yintercept = 0) +
    geom_point() +
    geom_text(
        aes(label = label1),
        colour = "black",
        size = 3,
        position = position_nudge(y = .05)
    ) +
    geom_text(
        aes(label = label2),
        colour = "black",
        size = 3,
        position = position_nudge(y = -.05)
    ) +
    scale_y_continuous(
        name = "",
        limits = c(-.15,.15), 
        expand = c(0,0),
        labels = NULL
    ) +
    scale_x_continuous(
        name = "",
        limits = c(0, 10),
        expand = c(0, 0),
        labels = NULL
    ) +
    scale_size_area(max_size = 15) +
    theme_minimal() +
    theme(
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none"
    )
```

```{r}
#| label: fig-international
#| fig-cap: "Proportionality of the allocation of ministries around the world. **Note**: Average marginal effect of seat share on ministry compensation depending on country. The closer to one, the more proportional the allocation of ministries."
#| fig-height: 3.37
#| fig-width: 6.1
#| fig-pos: t

mags <- left_join(
    avg_slopes(
        lm(
            min_share ~ cab_seat_share * country_name,
            df |>
                filter(
                    !is.na(cab_seat_share),
                    !is.na(min_share)
                ) |>
                filter(
                    coalition_type %in% c(2, 3, 5),
                    cabinet_party == 1
                ) |>
                group_by(cabinet_id) |>
                sample_n(size = n() - 1) |>
                ungroup() |>
                arrange(country_name_short, election_date)
        ),
        by = "country_name",
        variables = "cab_seat_share"
    ),
    df |>
        filter(
            !is.na(cab_seat_share),
            !is.na(min_share)
        ) |>
        filter(
            coalition_type %in% c(2, 3, 5),
            cabinet_party == 1
        ) |>
        group_by(cabinet_id) |>
        sample_n(size = n() - 1) |>
        ungroup() |>
        arrange(country_name_short, election_date) |>
        group_by(country_name) |>
        summarise(country_name_short = first(country_name_short))
        ) |>
        rename(sovereignt = country_name)

prop_conf <- mags |>
    filter(
        conf.high >= 1 & conf.low <=1
    ) |>
    nrow() |>
    spell_out()

cn_min1 <- mags$sovereignt[which.min(mags$estimate)]
cn_min2 <- mags$sovereignt[order(mags$estimate)[2]]

ggarrange(
    ggarrange(
        left_join(
            ne_countries(scale = 'medium', returnclass = 'sf'),
            mags,
            by = "sovereignt"
        ) |>
        ggplot(aes(fill = estimate)) +
            geom_sf(col = "grey30") +
            coord_sf(
                xlim = c(-25, 40),
                ylim = c(29, 73),
                expand = FALSE
            ) +
            scale_fill_gradient2(
                low = "navy",
                high = "firebrick",
                mid = "white",
                midpoint = 1,
                na.value = "grey80",
                limits = c(0, 2),
                name = "",
                guide = guide_colourbar(
                    frame.colour = 'black',
                    frame.linewidth = .5,
                    ticks.linewidth = 1,
                    ticks.colour = 'black'
                )
            ) +
            theme_bw() +
            theme(
                axis.text = element_blank(),
                axis.ticks = element_blank(),
                plot.margin = unit(c(2.8,2,0,2), 'mm'),
                legend.key.height = unit(2, 'mm'),
                legend.key.width = unit(15, 'mm'),
                legend.text = element_text(size=5),
                legend.title = element_text(size=8),
                panel.grid = element_blank()
            ),
        ggarrange(
            left_join(
                ne_countries(scale = "medium", returnclass = 'sf'),
                mags,
                by = "sovereignt"
            ) |>
            ggplot(aes(fill = estimate)) +
                geom_sf(col = "grey30") +
                coord_sf(
                    xlim = c(128, 146), 
                    ylim = c(28,46), 
                    expand = FALSE
                ) +
                scale_fill_gradient2(
                    low = "navy",
                    high = "firebrick",
                    mid = "white",
                    midpoint = 1,
                    na.value = "grey80",
                    limits = c(0, 2),
                    name = "",
                    guide = guide_colourbar(
                        frame.colour = 'black',
                        frame.linewidth = .5,
                        ticks.linewidth = 1,
                        ticks.colour = 'black'
                    )
                ) +
                scale_x_continuous(labels = NULL, breaks = NULL) +
                theme_bw() +
                theme(
                    axis.text = element_blank(),
                    axis.ticks = element_blank(),
                    plot.margin = unit(c(0,2,0,2), "mm"),
                    legend.key.height = unit(2, "mm"),
                    legend.key.width = unit(15, "mm"),
                    legend.text = element_text(size=5),
                    legend.title = element_text(size=8),
                    panel.grid = element_blank()
                ),
            left_join(
                ne_countries(scale = "medium", returnclass = 'sf'),
                mags,
                by = "sovereignt"
            ) |>
            ggplot(aes(fill = estimate)) +
                geom_sf(col = "grey30") +
                coord_sf(
                    xlim = c(112, 180), 
                    ylim = c(-48,-10), 
                    expand = FALSE
                ) +
                scale_fill_gradient2(
                    low = "navy",
                    high = "firebrick",
                    mid = "white",
                    midpoint = 1,
                    na.value = "grey80",
                    limits = c(0, 2),
                    name = "",
                    guide = guide_colourbar(
                        frame.colour = 'black',
                        frame.linewidth = .5,
                        ticks.linewidth = 1,
                        ticks.colour = 'black'
                    )
                ) +
                theme_bw() +
                theme(
                    axis.text = element_blank(),
                    axis.ticks = element_blank(),
                    plot.margin = unit(c(0,2,0,2), "mm"),
                    legend.key.height = unit(2, "mm"),
                    legend.key.width = unit(15, "mm"),
                    legend.text = element_text(size=5),
                    legend.title = element_text(size=8),
                    panel.grid = element_blank()
                ),
            ncol = 1,
            common.legend = F,
            legend = "none",
            heights = c(1.9, 1)
        ),
        ncol = 2,
        common.legend = T, 
        legend = 'bottom',
        widths = c(1.89, 1)
    ),
    mags |>
    ggplot(aes(fct_reorder(country_name_short, estimate), estimate)) +
        geom_hline(yintercept = 1, lwd = .3, linetype = "dashed", colour = "lightcoral") +
        geom_point(size = .8) +
        geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, lwd = .5) +
        coord_flip() +
        scale_x_discrete(name = "") +
        scale_y_continuous(
            name = "", 
            expand = c(0,0),
            limits = c(0, 2), 
            breaks = seq(0,2,.5),
        )  +
        theme_grey() +
        theme(
            axis.text.y = element_text(size = 5),
            axis.text.x = element_text(size = 8),
            plot.margin = unit(c(2,4,2,0), "mm"),
            panel.grid.major.x = element_blank(),
            panel.grid.minor.x = element_blank(),
            panel.border = element_rect(fill = NA, size = .8)
        ),
    ncol = 2,
    widths = c(2, 1)
)

country_n <- spell_out(nrow(mags))
```

```{r}
#| label: tbl-main-models
#| tbl-cap: Main selection models
#| results: "asis"
#| tbl-pos: t

if (knitr::is_html_output()){
    al <- "lcccc"
} else if (knitr::is_latex_output()){
    al <- "ldddd"
}

m1 <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r,
    data = df
)
m2 <- selection(
    cabinet_party ~ par_seat_share + left_right_r + country_name - 1, 
    dif_rou ~ cab_seat_share + left_right_r,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr - 1,
                make_matrix(df, 3, c("exp", "ps", "pr"))
            ),
            m1,
            c(1,2),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + rho,
                make_matrix(dfc, 4, c("exp", "cs", "pr", "rho"))
            ),
            m1,
            c(36:38),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr - 1,
                make_matrix(df, 3, c("exp", "ps", "lr"))
            ),
            m2,
            c(1,2),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + rho,
                make_matrix(dfc, 4, c("exp", "cs", "lr", "rho"))
            ),
            m2,
            c(36:38),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr" = "Relative radicalism",
        "lr" = "Absolute radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(m1, 13, 15),
        heck_gof(m2, 13, 15) |>
        select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = al,
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(10, hline_after = hl, extra_css = css) |>
row_spec(c(13:16), align = c("l", rep("c", 4)))


preds <- predictions(
        m1,
        newdata = datagrid(parlmean_dif_r = c(0, 4)),
        vcov = vcovCL(m1, cluster = df$country_name),
        type = "unconditional"
    )

midpos <- preds$estimate[1]
midse <- preds$std.error[1]
radpos <- preds$estimate[2]
radse <- preds$std.error[2]
pred_dif <- as.numeric(format(round(radpos - midpos, 2), nsmall = 2))
pred_lci <- format(round((radpos - midpos) - 1.96 * sqrt(radse^2 + midse^2), 2), nsmall = 2)
pred_uci <- format(round((radpos - midpos) + 1.96 * sqrt(radse^2 + midse^2), 2), nsmall = 2)
median_size <- round(
    df |> 
        group_by(cabinet_id) |> 
        summarise(ministries = max(cab_ministries)) |> 
        summarise(value = median(ministries, na.rm = TRUE)) |> 
        as.numeric()
)
dif_mini <- spell_out(round(sqrt((16 * (pred_dif / 100))^2)), max_value = 12)
```

```{r}
#| label: tbl-mechanism
#| tbl-cap: Selection models with both radicalism variants
#| results: "asis"
#| tbl-pos: t

if (knitr::is_html_output()){
    al <- "lcccc"
} else if (knitr::is_latex_output()){
    al <- "ldddd"
}

m3 <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + left_right_r + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + left_right_r,
    data = df
)
m4 <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r*left_right_r + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r*left_right_r,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + lr - 1,
                make_matrix(df, 4, c("exp", "ps", "pr", "lr"))
            ),
            m3,
            c(1:3),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + lr + rho,
                make_matrix(dfc, 5, c("exp", "cs", "pr", "lr", "rho"))
            ),
            m3,
            c(37:40),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + lr + prlr - 1,
                make_matrix(df, 5, c("exp", "ps", "pr", "lr", "prlr"))
            ),
            m4,
            c(1:3, 37),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + lr + prlr + rho,
                make_matrix(dfc, 6, c("exp", "cs", "pr", "lr", "prlr", "rho"))
            ),
            m4,
            c(38:42),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr" = "Relative radicalism",
        "lr" = "Absolute radicalism",
        "prlr" = "Absolute$\\times$relative radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(m3, 15, 17),
        heck_gof(m4, 15, 17) |>
        select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = al,
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(12, hline_after = hl, extra_css = css) |>
row_spec(c(15:18), align = c("l", rep("c", 4)))
```

```{r}
#| label: fig-ame
#| fig-cap: "Marginal Effects of Radicalism. **Note**: Each panel depicts its title's variable's predicted effect on a party's ministry compensation under different values of the other. 95 per cent confidence interval shown as ribbon. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1
#| fig-pos: t

rbind(
    slopes(
        m4,
        variables = "left_right_r",
        newdata = datagrid(
            parlmean_dif_r = seq(0, 5, .1)
        ),
        vcov = vcovCL(m4, cluster = df$country_name),
        type = "unconditional"
    ) |>
    rename(condition = parlmean_dif_r) |>
    select(-c(predicted:cab_seat_share)),
    slopes(
        m4,
        variables = "parlmean_dif_r",
        newdata = datagrid(
            left_right_r = seq(0, 5, .1)
        ),
        vcov = vcovCL(m4, cluster = df$country_name),
        type = "unconditional"
    ) |>
    rename(condition = left_right_r) |>
    select(-c(predicted:cab_seat_share))
) |>
mutate(term = ifelse(
        term == "left_right_r", 
        "Marginal effect of\nabsolute radicalism", 
        "Marginal effect of\nrelative radicalism"
    )
) |>
ggplot(aes(condition, estimate)) +
    facet_wrap(~term) +
    geom_hline(yintercept = 0, colour = "lightcoral") +
    geom_line() +
    geom_rug(
        data = dfc |>
            pivot_longer(
                c(parlmean_dif_r, left_right_r),
                names_to = "term",
                values_to = "condition"
            ) |>
            mutate(term = ifelse(
                term == "parlmean_dif_r", 
                "Marginal effect of\nabsolute radicalism", 
                "Marginal effect of\nrelative radicalism")
            ), 
        aes(x = condition),
        inherit.aes = FALSE,
        alpha = .2
    ) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3) +
    scale_y_continuous(name = "Marginal Effect") +
    scale_x_continuous(
        name = "Level of ...\n... relative radicalism                                         ... absolute radicalism",
        limits = c(0,5),
        expand = c(0,0)
    ) +
    theme(panel.spacing = unit(2, "lines"))
```

```{r}
#| label: three-step-preds

dfh <- df |>
    group_by(party_id) |>
    filter(
        !is.na(dif_rou_dif),
        previous_cabinet_id == lag(cabinet_id)
    ) |>
    ungroup()

dfhc <- dfh |> filter(cab_lag == 1)

sel <- selection(
        cab_lag ~ par_seat_share_l + parlmean_dif_rl + country_name - 1,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + parlmean_dif_rl + parlmean_dif_r_dif,
        data = dfh
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfh, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + parlmean_dif_rl + parlmean_dif_r_dif + dif_rou_l + imr,
    data = dfh2
)

ps1 <- predictions(
    sel,
    newdata = datagrid(
        par_seat_share_l = c(.2, .2),
        parlmean_dif_rl = c(0, 4)
    ),
    type = "response",
    vcov = vcovCL(sel, cluster = dfh$country_name)
)

cini <- paste0(round(ps1$estimate[1] * 100, 0), "%")
rini <- paste0(round(ps1$estimate[2] * 100, 0), "%")

ps2 <- predictions(
    sel,
    newdata = datagrid(
        par_seat_share_l = c(.2, .2, .2),
        par_seat_share_dif = c(0, 0, 0),
        parlmean_dif_rl = c(0, 4, 4),
        parlmean_dif_r_dif = c(0, 0, 1)
    ),
    type = "unconditional",
    vcov = vcovCL(sel, cluster = dfh$country_name)
) |>
mutate(estimate = pnorm(estimate))

ccon <- paste0(round(ps2$estimate[1] * 100, 0), "%")
r1con <- paste0(round(ps2$estimate[3] * 100, 0), "%")
r2con <- paste0(round(ps2$estimate[4] * 100, 0), "%")

ccomb <- paste0(round(round(ps1$estimate[1], 2) * round(ps2$estimate[1], 2) * 100, 0), "%")
r1comb <- paste0(round(round(ps1$estimate[2], 2) * round(ps2$estimate[3], 2) * 100, 0), "%")
r2comb <- paste0(round(round(ps1$estimate[2], 2) * round(ps2$estimate[4], 2) * 100, 0), "%")

psx <- predictions(
    m1, 
    newdata = datagrid(cab_seat_share = .2, parlmean_dif_r = c(0,4)),
    vcov = vcovCL(m1, df$country_name),
    type = "unconditional"
)


ps3 <- predictions(
    out,
    newdata = datagrid(
        cab_seat_share_l = c(1/3, 1/3, 1/3),
        cab_seat_share_dif = c(0, 0, 0),
        parlmean_dif_rl = c(0, 4, 4),
        parlmean_dif_r_dif = c(0, 0, 1),
        dif_rou_l = c(psx$estimate[1],psx$estimate[2],psx$estimate[2])
    ),
    vcov = vcovCL(out, cluster = dfh2$country_name)
)

cgam <- ifelse(round(ps3$estimate[1], 1) > 0, paste0("+", round(ps3$estimate[1], 1)), round(ps3$estimate[1], 1))
r1gam <- ifelse(round(ps3$estimate[6], 1) > 0, paste0("+", round(ps3$estimate[6], 1)), round(ps3$estimate[6], 1))
r2gam <- ifelse(round(ps3$estimate[8], 1) > 0, paste0("+", round(ps3$estimate[8], 1)), round(ps3$estimate[8], 1))
```

```{r}
#| label: fig-pred
#| fig-cap: "Predicted ministry compensation. **Note:** Thick line represents unconditional estimation. Shaded area shows 95% confidence interval of this estimation. Other lines represent estimations accounting for one of the considered moderators. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1
#| fig-pos: t

m1r <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE),
    data = df
)
m2r <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):auth_history + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):auth_history,
    data = df
)
m3r <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):cold_war + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):cold_war,
    data = df
)
m4r <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):year + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):year,
    data = df
)
m5r <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):time_per + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):time_per,
    data = df
)
m6r <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):eu_member + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):eu_member,
    data = df
)
m7r <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):seats_eu + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):seats_eu,
    data = df
)
m1a <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE),
    data = df
)
m2a <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):auth_history + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):auth_history,
    data = df
)
m3a <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):cold_war + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):cold_war,
    data = df
)
m4a <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):time_per + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):time_per,
    data = df
)
m5a <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):eu_member + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):eu_member,
    data = df
)
m6a <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):seats_eu + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):seats_eu,
    data = df
)

lr_pred <- predictions(
        m1r,
        newdata = datagrid(parlmean_dif = seq(-5,5, .1)),
        vcov = vcovCL(m1r, cluster = df$country_name),
        type = "unconditional"
    ) |>
    mutate(
        main = TRUE,
        rad = "Relative radicalism",
        int = "none"
    ) |>
    rename(radicalism = parlmean_dif) |>
    select(estimate, conf.low, conf.high, radicalism, main, rad, int) |>
    rbind(
        predictions(
            m2r,
            newdata = datagrid(
                parlmean_dif = rep(seq(-5, 5, .1), 3),
                auth_history = c(
                    rep("right", 101),
                    rep("left", 101),
                    rep("neither", 101)
                )
            ),
            vcov = vcovCL(m2r, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Relative radicalism",
            int = case_when(
                auth_history == "right" ~ "r_auth",
                auth_history == "left" ~ "l_auth",
                auth_history == "neither" ~ "n_auth",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = parlmean_dif) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m3r,
            newdata = datagrid(
                parlmean_dif = rep(seq(-5, 5, .1), 2),
                cold_war = c(
                    rep("pre", 101),
                    rep("post", 101)
                )
            ),
            vcov = vcovCL(m3r, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Relative radicalism",
            int = case_when(
                cold_war == "pre" ~ "prcw",
                cold_war == "post" ~ "pocw",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = parlmean_dif) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m4r,
            newdata = datagrid(
                parlmean_dif = rep(seq(-5, 5, .1), 3),
                year = c(
                    rep(1960, 101),
                    rep(1990, 101),
                    rep(2020, 101)
                )
            ),
            vcov = vcovCL(m4r, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Relative radicalism",
            int = factor(case_when(
                year == 1960 ~ "60y",
                year == 1990 ~ "90y",
                year == 2020 ~ "20y",
                TRUE ~ NA
            )),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = parlmean_dif) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m5r,
            newdata = datagrid(
                parlmean_dif = rep(seq(-5, 5, .1), 4),
                time_per = c(
                    rep("1960 and earlier", 101),
                    rep("1961-1980", 101),
                    rep("1981-2000", 101),
                    rep("2001-2020", 101)
                )
            ),
            vcov = vcovCL(m5r, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Relative radicalism",
            int = case_when(
                time_per == "1960 and earlier" ~ "60-",
                time_per == "1961-1980" ~ "6180",
                time_per == "1981-2000" ~ "8120",
                time_per == "2001-2020" ~ "0120",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = parlmean_dif) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m6r,
            newdata = datagrid(
                parlmean_dif = rep(seq(-5, 5, .1), 2),
                eu_member = c(
                    rep(TRUE, 101),
                    rep(FALSE, 101)
                )
            ),
            vcov = vcovCL(m6r, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Relative radicalism",
            int = case_when(
                eu_member == TRUE ~ "eum",
                eu_member == FALSE ~ "eux",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = parlmean_dif) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m7r,
            newdata = datagrid(
                parlmean_dif = rep(seq(-5, 5, .1), 2),
                seats_eu = c(
                    rep(TRUE, 101),
                    rep(FALSE, 101)
                )
            ),
            vcov = vcovCL(m7r, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Relative radicalism",
            int = case_when(
                seats_eu == TRUE ~ "ept",
                seats_eu == FALSE ~ "epf",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = parlmean_dif) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m1a,
            newdata = datagrid(left_right = seq(0,10, .1)),
            vcov = vcovCL(m1a, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = TRUE,
            rad = "Absolute radicalism",
            int = "none"
        ) |>
        rename(radicalism = left_right) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m2a,
            newdata = datagrid(
                left_right = rep(seq(0, 10, .1), 3),
                auth_history = c(
                    rep("right", 101),
                    rep("left", 101),
                    rep("neither", 101)
                )
            ),
            vcov = vcovCL(m2a, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Absolute radicalism",
            int = case_when(
                auth_history == "right" ~ "r_auth",
                auth_history == "left" ~ "l_auth",
                auth_history == "neither" ~ "n_auth",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = left_right) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m3a,
            newdata = datagrid(
                left_right = rep(seq(0, 10, .1), 2),
                cold_war = c(
                    rep("pre", 101),
                    rep("post", 101)
                )
            ),
            vcov = vcovCL(m3a, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Absolute radicalism",
            int = case_when(
                cold_war == "pre" ~ "prcw",
                cold_war == "post" ~ "pocw",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = left_right) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m4a,
            newdata = datagrid(
                left_right = rep(seq(0, 10, .1), 4),
                time_per = c(
                    rep("1960 and earlier", 101),
                    rep("1961-1980", 101),
                    rep("1981-2000", 101),
                    rep("2001-2020", 101)
                )
            ),
            vcov = vcovCL(m4a, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Absolute radicalism",
            int = case_when(
                time_per == "1960 and earlier" ~ "60-",
                time_per == "1961-1980" ~ "6180",
                time_per == "1981-2000" ~ "8120",
                time_per == "2001-2020" ~ "0120",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = left_right) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m5a,
            newdata = datagrid(
                left_right = rep(seq(0, 10, .1), 2),
                eu_member = c(
                    rep(TRUE, 101),
                    rep(FALSE, 101)
                )
            ),
            vcov = vcovCL(m5a, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Absolute radicalism",
            int = case_when(
                eu_member == TRUE ~ "eum",
                eu_member == FALSE ~ "eux",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = left_right) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    ) |>
    rbind(
        predictions(
            m6a,
            newdata = datagrid(
                left_right = rep(seq(0, 10, .1), 2),
                seats_eu = c(
                    rep(TRUE, 101),
                    rep(FALSE, 101)
                )
            ),
            vcov = vcovCL(m6a, cluster = df$country_name),
            type = "unconditional"
        ) |>
        mutate(
            main = FALSE,
            rad = "Absolute radicalism",
            int = case_when(
                seats_eu == TRUE ~ "ept",
                seats_eu == FALSE ~ "epf",
                TRUE ~ NA
            ),
            conf.low = NA,
            conf.high = NA
        ) |>
        rename(radicalism = left_right) |>
        select(estimate, conf.low, conf.high, radicalism, main, rad, int)
    )

ggplot(aes(x = radicalism, y = estimate, alpha = main, lwd = main, colour = int), data = lr_pred) +
    facet_wrap(~rad, scales = "free_x") +
    geom_line() +
    ggrepel::geom_text_repel(
        data = lr_pred |>
            mutate(keep = case_when(
                rad == "Absolute radicalism" & int == "60-" & radicalism == 8 ~ TRUE,
                rad == "Absolute radicalism" & int == "6180" & radicalism == 9 ~ TRUE,
                rad == "Absolute radicalism" & int == "prcw" & radicalism == 10 ~ TRUE,
                TRUE ~ FALSE
            )) |>
            filter(keep == TRUE) |>
            mutate(lab = c("pre Cold War end", "before 1961", "1961-1980")),
        aes(label = lab),
        nudge_x = -3.5
    ) +
    geom_rug(
        data = rbind(
            dfc |> 
                mutate(
                    radicalism = parlmean_dif, 
                    rad = "Relative radicalism"
                ) |>
                select(radicalism, rad),
            dfc |> 
                mutate(
                    radicalism = left_right, 
                    rad = "Absolute radicalism"
                ) |>
                select(radicalism, rad)
        ),
        aes(x = radicalism),
        inherit.aes = FALSE, alpha = .2
    ) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
    scale_x_continuous(
        name = "",
        expand = c(0,0)
    ) +
    scale_y_continuous(name = "Predicted ministry compensation") +
    scale_alpha_manual(values = c(.4, 1)) +
    scale_linewidth_manual(values = c(.5,2)) +
    scale_colour_manual(values = rep(lc, 17)) +
    theme(
        legend.position = "none",
        panel.spacing = unit(2, "lines")
    )
```

# Appendix

```{r}
#| label: appendix-setup

library(dplyr)
library(ggplot2)
library(tidyr)
library(ggridges)
library(latex2exp)
library(kableExtra)
library(ggdark)
library(modelsummary)
library(ggpubr)
library(forcats)
library(sampleSelection)
library(lmtest)
library(sandwich)
library(tibble)
library(marginaleffects)
library(lme4)
library(glmmTMB)
library(ggh4x)
library(ordinal)
library(broom)
library(stringr)

set.seed(2119)

df <- read.csv("df.csv")
dfc <- df |>
    filter(cabinet_party == 1)

if (knitr::is_html_output()){
    theme_set(dark_theme_bw())
} else if (knitr::is_latex_output()){
    theme_set(theme_bw())
}

ltop <- if (knitr::is_latex_output()){
    "scale_down"
} else {
    NULL
}

hl <- if (knitr::is_latex_output()){
    TRUE
} else {
    FALSE
}

css <- if (knitr::is_html_output()) {
    "border-bottom: 1px solid"
} else {
    NULL
}
```

```{r}
#| label: appendix-functions

lower_ci <- function(mean, se, n){
    lower_ci <- (mean - qt(1 - ((1 - 0.95) / 2), n - 1) * se)
}

upper_ci <- function(mean, se, n){
    upper_ci <- mean + qt(1 - ((1 - 0.95) / 2), n - 1) * se
}

se <- function(sd, n) {
    se <- sd / sqrt(n)
}

change_coefs <- function(pseudo, model, vars, type, cluster = df$country_name){

    cluster <- coeftest(model, vcovCL, cluster = cluster)

    if(type == "sel"){
        pseudo <- summary(pseudo)
        pseudo$coefficients[,1] <- summary(model)$estimate[vars, 1]
        pseudo$coefficients[,2] <- cluster[vars, 2]
        pseudo$coefficients[,4] <- cluster[vars, 4]
    } else if(type == "out"){
        pseudo <- summary(pseudo)
        pseudo$coefficients[,1] <- summary(model)$estimate[c(vars, nrow(summary(model)$estimate)), 1]
        pseudo$coefficients[,2] <- cluster[c(vars, nrow(summary(model)$estimate)), 2]
        pseudo$coefficients[,4] <- cluster[c(vars, nrow(summary(model)$estimate)), 4]
    }

    pseudo
}

make_matrix <- function(data, nvar, colnames){
    x  <- matrix(rnorm(nrow(data) * nvar), ncol = nvar) |>
        data.frame()
    colnames(x) <- colnames

    x
}

heck_gof <- function(model, first, last){
    gofs <- model |>
        glance() |>
        t() |>
        data.frame() |>
        rownames_to_column() |>
        rename(v1 = 2) |>
        filter(rowname %in% c("AIC", "logLik")) |>
        mutate(
            v2 = "",
            rowname = ifelse(rowname == "AIC", "AIC", "Log Likelihood")
        ) |>
        select(rowname, v2, v1) |>
        rbind(list("Country FEs", "Yes", NA))

    gofs
}

cn_exl_tbl <- function(countries){

    modelsummary(
        list(
            "Selection" = change_coefs(
                lm(
                    exp ~ ps + pr - 1,
                    make_matrix(
                        df |> filter(country_name_short != countries[1]),
                        3,
                        c("exp", "ps", "pr")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[1])
                ),
                c(1:2),
                type = "sel",
                cluster = df |> 
                    filter(country_name_short != countries[1]) |> 
                    select(country_name)
            ),
            "Outcome" = change_coefs(
                lm(
                    exp ~ cs + pr + rho,
                    make_matrix(
                        df |> filter(country_name_short != countries[1]),
                        4,
                        c("exp", "cs", "pr", "rho")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[1])
                ),
                c(35:37),
                type = "out",
                cluster = df |> 
                    filter(country_name_short != countries[1]) |> 
                    select(country_name)
            ),
            "Selection" = change_coefs(
                lm(
                    exp ~ ps + pr - 1,
                    make_matrix(
                        df |> filter(country_name_short != countries[2]),
                        3,
                        c("exp", "ps", "pr")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[2])
                ),
                c(1:2),
                type = "sel",
                cluster = df |> 
                    filter(country_name_short != countries[2]) |> 
                    select(country_name)
            ),
            "Outcome" = change_coefs(
                lm(
                    exp ~ cs + pr + rho,
                    make_matrix(
                        df |> filter(country_name_short != countries[2]),
                        4,
                        c("exp", "cs", "pr", "rho")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[2])
                ),
                c(35:37),
                type = "out",
                cluster = df |> 
                    filter(country_name_short != countries[2]) |> 
                    select(country_name)
            ),
            "Selection" = change_coefs(
                lm(
                    exp ~ ps + pr - 1,
                    make_matrix(
                        df |> filter(country_name_short != countries[3]),
                        3,
                        c("exp", "ps", "pr")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[3])
                ),
                c(1:2),
                type = "sel",
                cluster = df |> 
                    filter(country_name_short != countries[3]) |> 
                    select(country_name)
            ),
            "Outcome" = change_coefs(
                lm(
                    exp ~ cs + pr + rho,
                    make_matrix(
                        df |> filter(country_name_short != countries[3]),
                        4,
                        c("exp", "cs", "pr", "rho")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[3])
                ),
                c(35:37),
                type = "out",
                cluster = df |> 
                    filter(country_name_short != countries[3]) |> 
                    select(country_name)
            ),
            "Selection" = change_coefs(
                lm(
                    exp ~ ps + pr - 1,
                    make_matrix(
                        df |> filter(country_name_short != countries[4]),
                        3,
                        c("exp", "ps", "pr")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[4])
                ),
                c(1:2),
                type = "sel",
                cluster = df |> 
                    filter(country_name_short != countries[4]) |> 
                    select(country_name)
            ),
            "Outcome" = change_coefs(
                lm(
                    exp ~ cs + pr + rho,
                    make_matrix(
                        df |> filter(country_name_short != countries[4]),
                        4,
                        c("exp", "cs", "pr", "rho")
                    )
                ),
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[4])
                ),
                c(35:37),
                type = "out",
                cluster = df |> 
                    filter(country_name_short != countries[4]) |> 
                    select(country_name)
            )
        ),
        estimate = "{estimate}{stars}",
        stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
        coef_map = c(
            "ps" = "Parliamentary seat share",
            "cs" = "Cabinet seat share",
            "pr" = "Relative radicalism",
            "(Intercept)" = "Intercept",
            "rho" = "$\\rho$"
        ),
        gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
        add_rows = cbind(
            heck_gof(
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[1])
                ), 11, 13
            ),
            heck_gof(
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[2])
                ), 11, 13
            ) |>
            select(-rowname),
            heck_gof(
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[3])
                ), 11, 13
            ) |>
            select(-rowname),
            heck_gof(
                selection(
                    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                    dif_rou ~ cab_seat_share + parlmean_dif_r,
                    data = df |> filter(country_name_short != countries[4])
                ), 11, 13
            ) |>
            select(-rowname)
        ),
        notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
        align = if(knitr::is_latex_output()){
            "ldddddddd"
        } else if(knitr::is_html_output()){
            "lcccccccc"
        },
        output = "latex",
        escape = FALSE
    ) |>
    kable_styling(latex_options = ltop, full_width = FALSE) |>
    row_spec(0, bold = TRUE, align = "c") |>
    row_spec(8, hline_after = hl, extra_css = css) |>
    row_spec(c(11:14), align = c("l", rep("c", 8)))
}

imr_sel <- function(model){
    (dnorm(predict(model, part = "outcome", type = "unconditional")) / pnorm(predict(model, part = "outcome", type = "unconditional")))
}

change_coefs_lmout <- function(pseudo, model, vars, cluster = df$country_name){
    cluster <- coeftest(model, vcovCL, cluster = cluster)
    coef_imr <- coef(model)["imr"]
    se_imr <- cluster["imr", "Std. Error"]

    sigma <- summary(model)$sigma

    rho <- coef_imr / sigma
    se_rho <- se_imr / sigma

    t_rho <- rho / se_rho
    p_value <- 2 * pt(-abs(t_rho), df.residual(model))
    
    pseudo <- summary(pseudo)
    pseudo$coefficients[,1] <- append(
        cluster[c(vars, nrow(summary(model)$estimate)), 1], 
        rho
    )
    pseudo$coefficients[,2] <- append(
        cluster[c(vars, nrow(summary(model)$estimate)), 2],
        se_rho
    )
    pseudo$coefficients[,4] <- append(
        cluster[c(vars, nrow(summary(model)$estimate)), 4],
        p_value
    )

    pseudo
}

change_zero <- function(pseudo, model, type){
    pseudo <- summary(pseudo)
    model <- summary(model)
    if(type == "zero"){
        rownames(model$coefficients$zi) <- rownames(pseudo$coefficients)
        pseudo$coefficients <- model$coefficients$zi
    } else if (type == "conditional"){
        rownames(model$coefficients$cond) <- rownames(pseudo$coefficients)
        pseudo$coefficients <- model$coefficients$cond
    }

    pseudo
}


zero_gof <- function(model){
    gofs <- model |>
        summary()
    
    format(round(gofs$AICtab[1:3], 3), nsmall = 3) |>
        as.data.frame() |>
        rownames_to_column() |>
        rename(v1 = 2) |>
        mutate(
            v2 = "",
            rowname = ifelse(rowname == "logLik", "Log Likelihood", rowname)
        ) |>
        select(rowname, v2, v1) |>
        rbind(
            data.frame(
                rowname = c(
                    "Cabinets",
                    "Countries"
                ),
                v2 = gofs$ngrps$zi,
                v1 = c(
                    "",
                    gofs$ngrps$cond
                )
            )
        )
}

```

```{r}
#| label: sd-data

df_sd <- df |>
    group_by(country_name) |>
    mutate(
        reld1 = ifelse(
            parlmean_dif > mean(parlmean_dif, na.rm = TRUE) + sd(parlmean_dif, na.rm = TRUE) | parlmean_dif < mean(parlmean_dif, na.rm = TRUE) - sd(parlmean_dif, na.rm = TRUE), 
            TRUE, FALSE
        ),
        reld196 = ifelse(
            parlmean_dif > mean(parlmean_dif, na.rm = TRUE) + 1.96 * sd(parlmean_dif, na.rm = TRUE) | parlmean_dif < mean(parlmean_dif, na.rm = TRUE) - 1.96 - sd(parlmean_dif, na.rm = TRUE), 
            TRUE, FALSE
        ),
        absd1 = ifelse(
            left_right > mean(left_right, na.rm = TRUE) + sd(left_right, na.rm = TRUE) | left_right < mean(left_right, na.rm = TRUE) - sd(left_right, na.rm = TRUE), 
            TRUE, FALSE
        ),
        absd196 = ifelse(
            left_right > mean(left_right, na.rm = TRUE) + 1.96 * sd(left_right, na.rm = TRUE) | left_right < mean(left_right, na.rm = TRUE) - 1.96 * sd(left_right, na.rm = TRUE), 
            TRUE, FALSE
        )
    ) |>
    ungroup()
```

```{r}
#| label: fig-selbias
#| fig-cap: "Party Positions. **Note**: Party positions according to ParlGov and CHES data. Each panel shows the density of the positions of all parties compared to only government parties."
#| fig-height: 6
#| fig-width: 6.1

df |>
    select(left_right, lrgen, cabinet_party) |>
    pivot_longer(cols = c(left_right, lrgen)) |>
    filter(value >= 0) |>
    mutate(gov_value = ifelse(cabinet_party == 1, value, NA)) |>
    rename(data = name) |>
    pivot_longer(cols = c(value, gov_value)) |>
    filter(value >= 0) |>
    mutate(
        name = ifelse(name == "value", "All", "Government"),
        data = ifelse(data == "left_right", "ParlGov", "CHES")
    ) |>
    select(-cabinet_party) |>
    ggplot(aes(value, name, fill = stat(x))) +
        facet_wrap(~data, ncol = 1) +
        geom_density_ridges_gradient(
            bandwidth = .3,
            scale = 1.1,
            quantile_lines = TRUE,
            quantile_fun = function(x, ...)mean(x)
        ) +
        scale_fill_gradient2(
            low = "red3",
            high = "mediumblue",
            mid = "white",
            midpoint = 5,
            guide = "none"
        ) +
        scale_y_discrete(
            name = "Density",
            expand = c(0,0)
        ) +
        scale_x_continuous(
            name = "Left-right position",
            limits = c(0,10),
            expand = c(0,0),
            breaks = seq(0,10,1),
            labels = c("left", rep("", 4), "centre", rep("", 4), "right")
        ) +
        theme(
            plot.margin = unit(c(2, 6, 2, 2), "mm"),
            panel.spacing = unit(2, "lines")
        )
```

```{r}
#| label: fig-rads
#| fig-cap: "Party classification. **Note:** Party classification based on statistical criteria with standard deviations as red lines. PopuList coding as far left/right shown via shape, PopuList coding of classification as populist shown via colour. Cabinet status of parties indicated via transparency."
#| fig-height: 6
#| fig-width: 6.1

df |>
    group_by(country_name) |>
    mutate(
        rel_sd = (parlmean_dif - mean(parlmean_dif, na.rm = TRUE)) / sd(parlmean_dif, na.rm = TRUE),
        abs_sd = (left_right - mean(left_right, na.rm = TRUE)) / sd(left_right, na.rm = TRUE),
        pl = factor(case_when(
            pl_coded == FALSE ~ "not coded",
            populist == TRUE ~ "TRUE",
            populist == FALSE ~ "FALSE",
            TRUE ~ NA
        ), levels = c("FALSE", "TRUE", "not coded")),
        lr = factor(case_when(
            pl_coded == FALSE ~ "not coded",
            far_lr == TRUE ~ "TRUE",
            far_lr == FALSE ~ "FALSE",
            TRUE ~ NA
        ), levels = c("FALSE", "TRUE", "not coded"))
    ) |>
    ggplot(aes(rel_sd, abs_sd, alpha = as.logical(cabinet_party), colour = pl, shape = lr)) +
        geom_vline(
            xintercept = c(-1.96, -1, 0, 1, 1.96), 
            colour = c("firebrick", "lightcoral", "grey30", "lightcoral", "firebrick"), 
            alpha = .7
        ) +
        geom_hline(
            yintercept = c(-1.96, -1, 0, 1, 1.96), 
            colour = c("firebrick", "lightcoral", "grey30", "lightcoral", "firebrick"), 
            alpha = .7
        ) + 
        geom_point() +
        scale_x_continuous(
            name = TeX("$\\frac{RelativeRadicalism_{ik}-\\bar{RelativeRadicalism_{k}}}{\\sigma_{RelativeRadicalism_{k}}}$")
        ) +
        scale_y_continuous(
            name = TeX("$\\frac{LeftRight_{ik}-\\bar{LeftRight_{k}}}{\\sigma_{LeftRight_{k}}}$")
        ) +
        scale_alpha_discrete(name = "In government?") +
        scale_colour_manual(
            name = "PopuList: Populist?",
            values = c("red", "skyblue2", "grey30")
        ) +
        scale_shape_manual(
            name = "PopuList: Far left/right?",
            values = c(1, 18, 4)
        ) +
        guides(
            alpha=guide_legend(ncol=1, title.position = "top", title.hjust=0.5),
            shape=guide_legend(ncol=1, title.position = "top", title.hjust=0.5),
            colour=guide_legend(ncol=1, title.position = "top", title.hjust=0.5)
        ) +
        theme(
            legend.position = "bottom"
        )
```

```{r}
#| label: tbl-rads
#| tbl-cap: Radical parties in government

pinfo <- cbind(
    df_sd |>
        group_by(country_name) |>
        summarise(
            plc = max(pl_coded, na.rm = TRUE),
            across(
                c(reld1, reld196, absd1, absd196, populist, far_lr),
                ~ sum(., na.rm = TRUE)
            )
        ) |>
        mutate(
            across(
                c(populist, far_lr), 
                ~ ifelse(plc == FALSE, "not coded", as.character(.))
            )
        ) |>
        select(-plc),
    df_sd |>
        group_by(party_id) |>
        summarise(
            country_name = first(country_name),
            pl_coded = max(pl_coded, na.rm = TRUE),
            across(
                c(reld1, reld196, absd1, absd196, populist, far_lr),
                ~ max(., na.rm = TRUE)
            )
        ) |>
        group_by(country_name) |>
        summarise(
            plc = max(pl_coded, na.rm = TRUE),
            across(
                c(reld1, reld196, absd1, absd196, populist, far_lr),
                ~ sum(., na.rm = TRUE),
                .names = "{col}_un"
            )
        ) |>
        mutate(
            across(
                c(populist_un, far_lr_un), 
                ~ ifelse(plc == FALSE, "not coded", as.character(.))
            )
        ) |>
        select(-c(plc, country_name)),
    df_sd |>
        filter(cabinet_party == 1) |>
        group_by(country_name) |>
        summarise(
            plc = max(pl_coded, na.rm = TRUE),
            across(
                c(reld1, reld196, absd1, absd196, populist, far_lr),
                ~ sum(., na.rm = TRUE),
                .names = "{col}_gov"
            )
        ) |>
        mutate(
            across(
                c(populist_gov, far_lr_gov), 
                ~ ifelse(plc == FALSE, "not coded", as.character(.))
            )
        ) |>
        select(-c(plc, country_name)),
    df_sd |>
        filter(cabinet_party == 1) |>
        group_by(party_id) |>
        summarise(
            country_name = first(country_name),
            pl_coded = max(pl_coded, na.rm = TRUE),
            across(
                c(reld1, reld196, absd1, absd196, populist, far_lr),
                ~ max(., na.rm = TRUE)
            )
        ) |>
        group_by(country_name) |>
        summarise(
            plc = max(pl_coded, na.rm = TRUE),
            across(
                c(reld1, reld196, absd1, absd196, populist, far_lr),
                ~ sum(., na.rm = TRUE),
                .names = "{col}_gov_un"
            )
        ) |>
        mutate(
            across(
                c(populist_gov_un, far_lr_gov_un), 
                ~ ifelse(plc == FALSE, "not coded", as.character(.))
            )
        ) |>
        select(-c(plc, country_name))
)

pinfo_t <- pinfo |>
    summarise(across(everything(), ~ sum(as.numeric(.), na.rm = TRUE))) |>
    mutate(country_name = "Total") |>
    mutate(across(c(populist, far_lr, populist_un, populist_gov, populist_gov_un, far_lr_un, far_lr_gov, far_lr_gov_un), ~ as.character(.)))

bind_rows(pinfo, pinfo_t) |>
mutate(
    rd1 = paste0(reld1, "(", reld1_un, ") [", reld1_gov, "(", reld1_gov_un, ")]"),
    rd196 = paste0(reld196, "(", reld196_un, ") [", reld196_gov, "(", reld196_gov_un, ")]"),
    ad1 = paste0(absd1, "(", absd1_un, ") [", absd1_gov, "(", absd1_gov_un, ")]"),
    ad196 = paste0(absd196, "(", absd196_un, ") [", absd196_gov, "(", absd196_gov_un, ")]"),
    pl = ifelse(populist == "not coded", "not coded", paste0(populist, "(", populist_un, ") [", populist_gov, "(", populist_gov_un, ")]")),
    flr = ifelse(far_lr == "not coded", "not coded", paste0(far_lr, "(", far_lr_un, ") [", far_lr_gov, "(", far_lr_gov_un, ")]"))
) |>
select(country_name, rd1:flr) |>
kbl(
    col.names = c(
        "Country",
        "$x_{ik}\\notin(\\bar{x_{k}}-\\sigma_{x_{k}}, \\bar{x_{k}}+\\sigma_{x_{k}})$",
        "$x_{ik}\\notin(\\bar{x_{k}}-1.96\\sigma_{x_{k}}, \\bar{x_{k}}+1.96\\sigma_{x_{k}})$",
        "$x_{ik}\\notin(\\bar{x_{k}}-\\sigma_{x_{k}}, \\bar{x_{k}}+\\sigma_{x_{k}})$",
        "$x_{ik}\\notin(\\bar{x_{k}}-1.96\\sigma_{x_{k}}, \\bar{x_{k}}+1.96\\sigma_{x_{k}})$",
        "Populist",
        "Far left/right"
    ),
    align = "lcccccc",
    output = "latex",
    escape = FALSE,
    booktabs = TRUE,
    linesep = ""
) |>
kable_styling(latex_options = ltop) |>
add_header_above(c("", "Relative Radicalism" = 2, "Absolute Radicalism" = 2, "PopuList" = 2)) |>
row_spec(34, bold = TRUE) |>
row_spec(33, hline_after = hl, extra_css = css) |>
footnote(general = "Observations(Unique Parties) [Observations in Government(Unique Parties in Government)]", footnote_as_chunk = T, general_title = "Cell Format: ")
```

```{r}
#| label: deviations

lib_dev <- df |>
    filter(cabinet_party == 1 & family_name != "to be coded") |>
    group_by(family_name) |>
    summarise(
        left_right = mean(left_right, na.rm = TRUE),
        mean = mean(dif_rou, na.rm = TRUE),
        sd = sd(dif_rou, na.rm = TRUE),
        sum = sum(!is.na(dif_rou))
    ) |>
    mutate(
        se = se(sd, sum),
        lci = lower_ci(mean, se, sum),
        uci = upper_ci(mean, se, sum),
        family_name = reorder(family_name, left_right)
    ) |>
    filter(family_name == "Liberal") |>
    select(mean) |>
    as.numeric() |>
    round(digits = 2) |>
    format(nsmall = 2)

min_cab_lib <- round(1 / (as.numeric(lib_dev) / 100))
```

```{r}
#| label: fig-descr
#| fig-cap: "Ministry compensation tendencies of the party families. **Note**: Party family-wise average ministry compensation. 95 per cent confidence intervals shown."
#| fig-height: 4.96
#| fig-width: 6.1

ggarrange(
    df |>
        group_by(party_id) |>
        summarise(
            left_right = first(left_right),
            family_name = first(family_name)
        ) |>
        mutate(family_name = reorder(family_name, left_right, FUN = mean)) |>
        filter(family_name != "to be coded") |>
        ggplot(aes(left_right, fct_rev(family_name), fill = stat(x))) +
            geom_density_ridges_gradient(
                bandwidth = .3,
                scale = 1.5,
                colour = "black",
                quantile_lines = TRUE,
                quantile_fun = function(x, ...)mean(x)
            ) +
            scale_fill_gradient2(
                low = "red3",
                high = "mediumblue",
                mid = "white",
                midpoint = 5,
                guide = "none"
            ) +
            scale_y_discrete(
                name = "Density",
                expand = c(0,0)
            ) +
            scale_x_continuous(
                name = "Position",
                limits = c(0,10),
                expand = c(0,0),
                breaks = seq(0,10,1),
                labels = c("left", rep("", 4), "centre", rep("", 4), "right")
            ) +
            theme(plot.margin = unit(c(2,6,4,2), "mm")),
    df |>
        filter(cabinet_party == 1 & family_name != "to be coded") |>
        group_by(family_name) |>
        summarise(
            left_right = mean(left_right, na.rm = TRUE),
            mean = mean(dif_rou, na.rm = TRUE),
            sd = sd(dif_rou, na.rm = TRUE),
            sum = sum(!is.na(dif_rou))
        ) |>
        mutate(
            se = se(sd, sum),
            lci = lower_ci(mean, se, sum),
            uci = upper_ci(mean, se, sum),
            family_name = reorder(family_name, left_right)
        ) |>
        ggplot(aes(mean, fct_rev(family_name), colour = left_right)) +
            geom_vline(
                xintercept = 0, 
                col = "lightcoral", 
                linetype = "longdash"
            ) +
            geom_point(size = 3) +
            geom_errorbarh(
                aes(xmin = lci, xmax = uci),
                height = 0, 
                linewidth = .8
            ) +
            geom_text(
                aes(label = paste(fct_rev(family_name), ":", sep = "")),
                x = -3.8,
                hjust = 0,
                vjust = -.7,
                col = "grey30",
                size = 3
            ) +
            scale_colour_gradient2(
                low = "red3",
                high = "mediumblue",
                mid = "grey80",
                midpoint = 5,
                guide = "none"
            ) +
            scale_y_discrete(
                name = "",
                expand = expansion(add = c(.15, .4))
            ) +
            scale_x_continuous(
                name = "Deviation",
                expand = c(0,0),
                limits = c(-4, 4)
            ) +
            theme(
                plot.margin = unit(c(14,4,2,2), 'mm'),
                axis.title = element_text(size = 10),
                axis.text.y = element_blank(),
                panel.grid.minor = element_blank()
            ),
    ncol = 2
)
```

```{r}
#| label: fig-case-sel
#| fig-caption: "Selection of mini case studies in the article. **Note**: Selected cases shown as crosses."
#| fig-height: 6
#| fig-width: 6.1

csm <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + left_right_r + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + left_right_r,
    data = df
)

df |>
    mutate(
        preds = predict(csm),
        sel = case_when(
            cabinet_id == 672 & party_id == 373 ~ TRUE,
            cabinet_id == 519 & party_id == 521 ~ TRUE,
            cabinet_id == 1079 & party_id == 2268 ~ TRUE,
            cabinet_id == 1510 & party_id == 2645 ~ TRUE,
            TRUE ~ FALSE
        ),
        residual = residuals(csm)
    ) |>
    filter(cabinet_party == 1) |>
    ggplot(aes(left_right_r, residual, colour = parlmean_dif_r, size = sel, shape = sel, label = party_name)) +
        geom_point(stroke = 3) +
        scale_x_continuous(
            name = "Absolute radicalism",
            limits = c(0, 5),
            expand = c(0,0)
        ) +
        scale_y_continuous(name = "Residual") +
        scale_colour_gradient(
            name = "Relative radicalism",
            limits = c(0,5),
            low = "yellow",
            high = "red"
        ) +
        scale_shape_manual(
            name = "",
            values = c("TRUE" = 4, "FALSE" = 16),
            labels = c("TRUE" = "Selected", "FALSE" = "Not selected")
        ) +
        scale_size_manual(
            name = "",
            values = c("TRUE" = 10, "FALSE" = 2),
            labels = c("TRUE" = "Selected", "FALSE" = "Not selected")
        ) +
        theme_bw() +
        theme(legend.position = "bottom")
```

```{r}
#| label: ches_data

df_ches <- df |>
    filter(
        !is.na(lrgen_r),
        !is.na(lrecon_r),
        !is.na(galtan_r)
    )

df_ches_c <- df_ches |>
    filter(cabinet_party == 1)
```

```{r}
#| label: tbl-ches
#| tbl-cap: Heckman models with CHES data
#| results: "asis"

m1ches <- selection(
    cabinet_party ~ par_seat_share + lrgen_r + country_name - 1, 
    dif_rou ~ cab_seat_share + lrgen_r,
    data = df_ches
)
m2ches <- selection(
    cabinet_party ~ par_seat_share + lrecon_r + country_name - 1, 
    dif_rou ~ cab_seat_share + lrecon_r,
    data = df_ches
)
m3ches <- selection(
    cabinet_party ~ par_seat_share + galtan_r + country_name - 1, 
    dif_rou ~ cab_seat_share + galtan_r,
    data = df_ches
)
m4ches <- selection(
    cabinet_party ~ par_seat_share + lrecon_r + galtan_r + country_name - 1, 
    dif_rou ~ cab_seat_share + lrecon_r + galtan_r,
    data = df_ches
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr - 1,
                make_matrix(df_ches, 3, c("exp", "ps", "lr"))
            ),
            m1ches,
            c(1,2),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + rho,
                make_matrix(df_ches_c, 4, c("exp", "cs", "lr", "rho"))
            ),
            m1ches,
            c(28:30),
            type = "out",
            cluster = df_ches$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec - 1,
                make_matrix(df_ches, 3, c("exp", "ps", "ec"))
            ),
            m2ches,
            c(1,2),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + rho,
                make_matrix(df_ches_c, 4, c("exp", "cs", "ec", "rho"))
            ),
            m2ches,
            c(28:30),
            type = "out",
            cluster = df_ches$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + gt - 1,
                make_matrix(df_ches, 3, c("exp", "ps", "gt"))
            ),
            m3ches,
            c(1,2),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + gt + rho,
                make_matrix(df_ches_c, 4, c("exp", "cs", "gt", "rho"))
            ),
            m3ches,
            c(28:30),
            type = "out",
            cluster = df_ches$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + gt - 1,
                make_matrix(df_ches, 4, c("exp", "ps", "ec", "gt"))
            ),
            m4ches,
            c(1:3),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + gt + rho,
                make_matrix(df_ches_c, 5, c("exp", "cs", "ec", "gt", "rho"))
            ),
            m4ches,
            c(29:32),
            type = "out",
            cluster = df_ches$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "lr" = "Left-right radicalism",
        "ec" = "Economic radicalism",
        "gt" = "Cultural radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(m1ches, 17, 19),
        heck_gof(m2ches, 17, 19) |> select(-rowname),
        heck_gof(m3ches, 17, 19) |> select(-rowname),
        heck_gof(m4ches, 17, 19) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(12, hline_after = hl, extra_css = css) |>
row_spec(c(15:18), align = c("l", rep("c", 8)))
```

```{r}
#| label: tbl-ches-side
#| tbl-cap: Heckman models with CHES data accounting for side of radicalism
#| results: "asis"

m1aches <- selection(
    cabinet_party ~ par_seat_share + poly(lrgen, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(lrgen, 2, raw = TRUE),
    data = df_ches
)
m2aches <- selection(
    cabinet_party ~ par_seat_share + poly(lrecon, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(lrecon, 2, raw = TRUE),
    data = df_ches
)
m3aches <- selection(
    cabinet_party ~ par_seat_share + poly(galtan, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(galtan, 2, raw = TRUE),
    data = df_ches
)
m4aches <- selection(
    cabinet_party ~ par_seat_share + poly(lrecon, 2, raw = TRUE) + poly(galtan, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(lrecon, 2, raw = TRUE) + poly(galtan, 2, raw = TRUE),
    data = df_ches
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr + lr2 - 1,
                make_matrix(df_ches, 4, c("exp", "ps", "lr", "lr2"))
            ),
            m1aches,
            c(1:3),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + lr2 + rho,
                make_matrix(df_ches_c, 5, c("exp", "cs", "lr", "lr2", "rho"))
            ),
            m1aches,
            c(29:32),
            type = "out",
            cluster = df_ches$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + ec2 - 1,
                make_matrix(df_ches, 4, c("exp", "ps", "ec", "ec2"))
            ),
            m2aches,
            c(1:3),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + ec2 + rho,
                make_matrix(df_ches_c, 5, c("exp", "cs", "ec", "ec2", "rho"))
            ),
            m2aches,
            c(29:32),
            type = "out",
            cluster = df_ches$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + gt + gt2 - 1,
                make_matrix(df_ches, 4, c("exp", "ps", "gt", "gt2"))
            ),
            m3aches,
            c(1:3),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + gt + gt2 + rho,
                make_matrix(df_ches_c, 5, c("exp", "cs", "gt", "gt2", "rho"))
            ),
            m3aches,
            c(29:32),
            type = "out",
            cluster = df_ches$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + ec2 + gt + gt2 - 1,
                make_matrix(df_ches, 6, c("exp", "ps", "ec", "ec2", "gt", "gt2"))
            ),
            m4aches,
            c(1:5),
            type = "sel",
            cluster = df_ches$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + ec2 + gt + gt2 + rho,
                make_matrix(df_ches_c, 7, c("exp", "cs", "ec", "ec2", "gt", "gt2", "rho"))
            ),
            m4aches,
            c(31:36),
            type = "out",
            cluster = df_ches$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "lr" = "Left-right position",
        "lr2" = "Left-right position$^{2}$",
        "ec" = "Economic position",
        "ec2" = "Economic position$^{2}$",
        "gt" = "Cultural position",
        "gt2" = "Cultural position$^{2}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(m1aches, 23, 25),
        heck_gof(m2aches, 23, 25) |> select(-rowname),
        heck_gof(m3aches, 23, 25) |> select(-rowname),
        heck_gof(m4aches, 23, 25) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(18, hline_after = hl, extra_css = css) |>
row_spec(c(21:24), align = c("l", rep("c", 8)))
```

```{r}
#| label: cmp-data

df_cmp <- df |>
    filter(
        !is.na(rile_r),
        !is.na(market_nat_r),
        !is.na(culture_nat_r)
    )

df_cmp_c <- df_cmp |>
    filter(cabinet_party == 1)
```

```{r}
#| label: tbl-cmp
#| tbl-cap: Heckman models with CMP data
#| results: "asis"

m1cmp <- selection(
    cabinet_party ~ par_seat_share + rile_r + country_name - 1, 
    dif_rou ~ cab_seat_share + rile_r,
    data = df_cmp
)
m2cmp <- selection(
    cabinet_party ~ par_seat_share + market_nat_r + country_name - 1, 
    dif_rou ~ cab_seat_share + market_nat_r,
    data = df_cmp
)
m3cmp <- selection(
    cabinet_party ~ par_seat_share + culture_nat_r + country_name - 1, 
    dif_rou ~ cab_seat_share + culture_nat_r,
    data = df_cmp
)
m4cmp <- selection(
    cabinet_party ~ par_seat_share + market_nat_r + culture_nat_r + country_name - 1, 
    dif_rou ~ cab_seat_share + market_nat_r + culture_nat_r,
    data = df_cmp
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr - 1,
                make_matrix(df_cmp, 3, c("exp", "ps", "lr"))
            ),
            m1cmp,
            c(1,2),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + rho,
                make_matrix(df_cmp_c, 4, c("exp", "cs", "lr", "rho"))
            ),
            m1cmp,
            c(36:38),
            type = "out",
            cluster = df_cmp$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec - 1,
                make_matrix(df_cmp, 3, c("exp", "ps", "ec"))
            ),
            m2cmp,
            c(1,2),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + rho,
                make_matrix(df_cmp_c, 4, c("exp", "cs", "ec", "rho"))
            ),
            m2cmp,
            c(36:38),
            type = "out",
            cluster = df_cmp$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + gt - 1,
                make_matrix(df_cmp, 3, c("exp", "ps", "gt"))
            ),
            m3cmp,
            c(1,2),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + gt + rho,
                make_matrix(df_cmp_c, 4, c("exp", "cs", "gt", "rho"))
            ),
            m3cmp,
            c(36:38),
            type = "out",
            cluster = df_cmp$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + gt - 1,
                make_matrix(df_cmp, 4, c("exp", "ps", "ec", "gt"))
            ),
            m4cmp,
            c(1:3),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + gt + rho,
                make_matrix(df_cmp_c, 5, c("exp", "cs", "ec", "gt", "rho"))
            ),
            m4cmp,
            c(37:40),
            type = "out",
            cluster = df_cmp$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "lr" = "Left-right radicalism",
        "ec" = "Economic radicalism",
        "gt" = "Cultural radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(m1cmp, 17, 19),
        heck_gof(m2cmp, 17, 19) |> select(-rowname),
        heck_gof(m3cmp, 17, 19) |> select(-rowname),
        heck_gof(m4cmp, 17, 19) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(12, hline_after = hl, extra_css = css)  |>
row_spec(c(15:18), align = c("l", rep("c", 8)))
```

```{r}
#| label: tbl-cmp-side
#| tbl-cap: Heckman models with CMP data accounting for side of radicalism
#| results: "asis"

m1acmp <- selection(
    cabinet_party ~ par_seat_share + poly(rile, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(rile, 2, raw = TRUE),
    data = df_cmp
)
m2acmp <- selection(
    cabinet_party ~ par_seat_share + poly(market_nat, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(market_nat, 2, raw = TRUE),
    data = df_cmp
)
m3acmp <- selection(
    cabinet_party ~ par_seat_share + poly(culture_nat, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(culture_nat, 2, raw = TRUE),
    data = df_cmp
)
m4acmp <- selection(
    cabinet_party ~ par_seat_share + poly(market_nat, 2, raw = TRUE) + poly(culture_nat, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(market_nat, 2, raw = TRUE) + poly(culture_nat, 2, raw = TRUE),
    data = df_cmp
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr + lr2 - 1,
                make_matrix(df_cmp, 4, c("exp", "ps", "lr", "lr2"))
            ),
            m1acmp,
            c(1:3),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + lr2 + rho,
                make_matrix(df_cmp_c, 5, c("exp", "cs", "lr", "lr2", "rho"))
            ),
            m1acmp,
            c(37:40),
            type = "out",
            cluster = df_cmp$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + ec2 - 1,
                make_matrix(df_cmp, 4, c("exp", "ps", "ec", "ec2"))
            ),
            m2acmp,
            c(1:3),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + ec2 + rho,
                make_matrix(df_cmp_c, 5, c("exp", "cs", "ec", "ec2", "rho"))
            ),
            m2acmp,
            c(37:40),
            type = "out",
            cluster = df_cmp$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + gt + gt2 - 1,
                make_matrix(df_cmp, 4, c("exp", "ps", "gt", "gt2"))
            ),
            m3acmp,
            c(1:3),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + gt + gt2 + rho,
                make_matrix(df_cmp_c, 5, c("exp", "cs", "gt", "gt2", "rho"))
            ),
            m3acmp,
            c(37:40),
            type = "out",
            cluster = df_cmp$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + ec2 + gt + gt2 - 1,
                make_matrix(df_cmp, 6, c("exp", "ps", "ec", "ec2", "gt", "gt2"))
            ),
            m4acmp,
            c(1:5),
            type = "sel",
            cluster = df_cmp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + ec2 + gt + gt2 + rho,
                make_matrix(df_cmp_c, 7, c("exp", "cs", "ec", "ec2", "gt", "gt2", "rho"))
            ),
            m4acmp,
            c(39:44),
            type = "out",
            cluster = df_cmp$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "lr" = "Left-right position",
        "lr2" = "Left-right position$^{2}$",
        "ec" = "Economic position",
        "ec2" = "Economic position$^{2}$",
        "gt" = "Cultural position",
        "gt2" = "Cultural position$^{2}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(m1acmp, 23, 25),
        heck_gof(m2acmp, 23, 25) |> select(-rowname),
        heck_gof(m3acmp, 23, 25) |> select(-rowname),
        heck_gof(m4acmp, 23, 25) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(18, hline_after = hl, extra_css = css) |>
row_spec(c(21:24), align = c("l", rep("c", 8)))
```

```{r}
#| label: data-sq

dfsq <- df |>
    filter(
        !is.na(sqdif_dy),
        !is.na(sqdif_dy10),
        !is.na(sqdif_dy70),
        !is.na(sqdif_eq),
        !is.na(lastcab_dif)
    )

dfsq_c <- df |>
    filter(cabinet_party == 1)
```

```{r}
#| label: tbl-sq
#| tbl-cap: Heckman models with status quo distance
#| results: "asis"

msq1 <- selection(
    cabinet_party ~ par_seat_share + poly(sqdif_dy, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(sqdif_dy, 2, raw = TRUE),
    data = dfsq
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + sq + sq2 - 1,
                make_matrix(dfsq, 4, c("exp", "ps", "sq", "sq2"))
            ),
            msq1,
            c(1:3),
            type = "sel",
            cluster = dfsq$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + sq + sq2 + rho,
                make_matrix(dfsq_c, 5, c("exp", "cs", "sq", "sq2", "rho"))
            ),
            msq1,
            c(37:40),
            type = "out",
            cluster = dfsq$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "sq" = "Status quo distance",
        "sq2" = "Status quo distance²",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = heck_gof(msq1, 13, 15),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldd"
    } else if(knitr::is_html_output()){
        "lcc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = TRUE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(10, hline_after = hl, extra_css = css) |>
row_spec(c(13:16), align = c("l", rep("c", 2)))

```

```{r}
#| label: tbl-sq-alt
#| tbl-cap: Heckman models with alternative status quo distance
#| results: "asis"

msq2 <- selection(
    cabinet_party ~ par_seat_share + poly(sqdif_dy10, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(sqdif_dy10, 2, raw = TRUE),
    data = dfsq
)
msq3 <- selection(
    cabinet_party ~ par_seat_share + poly(sqdif_dy70, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(sqdif_dy70, 2, raw = TRUE),
    data = dfsq
)
msq4 <- selection(
    cabinet_party ~ par_seat_share + poly(sqdif_eq, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(sqdif_eq, 2, raw = TRUE),
    data = dfsq
)
msq5 <- selection(
    cabinet_party ~ par_seat_share + poly(lastcab_dif, 2, raw = TRUE) + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(lastcab_dif, 2, raw = TRUE),
    data = dfsq
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + sq10 + sq102 - 1,
                make_matrix(dfsq, 4, c("exp", "ps", "sq10", "sq102"))
            ),
            msq2,
            c(1:3),
            type = "sel",
            cluster = dfsq$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + sq10 + sq102 + rho,
                make_matrix(dfsq_c, 5, c("exp", "cs", "sq10", "sq102", "rho"))
            ),
            msq2,
            c(37:40),
            type = "out",
            cluster = dfsq$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + sq70 + sq702 - 1,
                make_matrix(dfsq, 4, c("exp", "ps", "sq70", "sq702"))
            ),
            msq3,
            c(1:3),
            type = "sel",
            cluster = dfsq$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + sq70 + sq702 + rho,
                make_matrix(dfsq_c, 5, c("exp", "cs", "sq70", "sq702", "rho"))
            ),
            msq3,
            c(37:40),
            type = "out",
            cluster = dfsq$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + sqe + sqe2 - 1,
                make_matrix(dfsq, 4, c("exp", "ps", "sqe", "sqe2"))
            ),
            msq4,
            c(1:3),
            type = "sel",
            cluster = dfsq$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + sqe + sqe2 + rho,
                make_matrix(dfsq_c, 5, c("exp", "cs", "sqe", "sqe2", "rho"))
            ),
            msq4,
            c(37:40),
            type = "out",
            cluster = dfsq$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lc + lc2 - 1,
                make_matrix(dfsq, 4, c("exp", "ps", "lc", "lc2"))
            ),
            msq5,
            c(1:3),
            type = "sel",
            cluster = dfsq$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lc + lc2 + rho,
                make_matrix(dfsq_c, 5, c("exp", "cs", "lc", "lc2", "rho"))
            ),
            msq5,
            c(37:40),
            type = "out",
            cluster = dfsq$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "sq10" = "SQD (10$\\%$)",
        "sq102" = "SQD (10$\\%$)$^{2}$",
        "sq70" = "SQD (70$\\%$)",
        "sq702" = "SQD (70$\\%$)$^{2}$",
        "sqe" = "SQD (0$\\%$)",
        "sqe2" = "SQD (0$\\%$)$^{2}$",
        "lc" = "Distance to last cabinet",
        "lc2" = "Distance to last cabinet$^{2}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(msq2, 25, 27),
        heck_gof(msq3, 25, 27) |> select(-rowname),
        heck_gof(msq4, 25, 27) |> select(-rowname),
        heck_gof(msq5, 25, 27) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(22, hline_after = hl, extra_css = css) |>
row_spec(c(25:28), align = c("l", rep("c", 8)))
```

```{r}
dfpl <- df |>
    filter(pl_coded == TRUE)

dfplc <- dfpl |>
    filter(cabinet_party == TRUE)

df_fam <- df |>
    filter(family_name != "to be coded")

df_fam_c <- df_fam |>
    filter(cabinet_party == 1)
```

```{r}
#| label: tbl-pl
#| tbl-cap: Heckman models with codings from PopuList
#| results: "asis"

mpl1 <- selection(
    cabinet_party ~ par_seat_share + country_name + populist - 1, 
    dif_rou ~ cab_seat_share + populist,
    data = dfpl
)
mpl2 <- selection(
    cabinet_party ~ par_seat_share + country_name + far_lr - 1, 
    dif_rou ~ cab_seat_share + far_lr,
    data = dfpl
)
mpl3 <- selection(
    cabinet_party ~ par_seat_share + country_name + far_left + far_right - 1, 
    dif_rou ~ cab_seat_share + far_left + far_right,
    data = dfpl
)
mpl4 <- selection(
    cabinet_party ~ par_seat_share + country_name + populist + far_left + far_right - 1, 
    dif_rou ~ cab_seat_share + populist + far_left + far_right,
    data = dfpl
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + po - 1,
                make_matrix(dfpl, 3, c("exp", "ps", "po"))
            ),
            mpl1,
            c(1,31),
            type = "sel",
            cluster = dfpl$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + po + rho,
                make_matrix(dfplc, 4, c("exp", "cs", "po", "rho"))
            ),
            mpl1,
            c(32:34),
            type = "out",
            cluster = dfpl$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + flr - 1,
                make_matrix(dfpl, 3, c("exp", "ps", "flr"))
            ),
            mpl2,
            c(1,31),
            type = "sel",
            cluster = dfpl$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + flr + rho,
                make_matrix(dfplc, 4, c("exp", "cs", "flr", "rho"))
            ),
            mpl2,
            c(32:34),
            type = "out",
            cluster = dfpl$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + fl + fr - 1,
                make_matrix(dfpl, 4, c("exp", "ps", "fl", "fr"))
            ),
            mpl3,
            c(1,31,32),
            type = "sel",
            cluster = dfpl$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + fl + fr + rho,
                make_matrix(dfplc, 5, c("exp", "cs", "fl", "fr", "rho"))
            ),
            mpl3,
            c(33:36),
            type = "out",
            cluster = dfpl$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + po + fl + fr - 1,
                make_matrix(dfpl, 5, c("exp", "ps", "po", "fl", "fr"))
            ),
            mpl4,
            c(1,31:33),
            type = "sel",
            cluster = dfpl$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + po + fl + fr + rho,
                make_matrix(dfplc, 6, c("exp", "cs", "po", "fl", "fr", "rho"))
            ),
            mpl4,
            c(34:38),
            type = "out",
            cluster = dfpl$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "po" = "Populist",
        "flr" = "Far left/right",
        "fl" = "Far left",
        "fr" = "Far right",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mpl1, 19, 21),
        heck_gof(mpl2, 19, 21) |> select(-rowname),
        heck_gof(mpl3, 19, 21) |> select(-rowname),
        heck_gof(mpl4, 19, 21) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(14, hline_after = hl, extra_css = css) |>
row_spec(c(17:20), align = c("l", rep("c", 8)))
```

```{r}
#| label: tbl-sd
#| tbl-cap: Heckman models with standard deviation based radicalism coding
#| results: "asis"

msd1 <- selection(
    cabinet_party ~ par_seat_share + country_name + reld1 - 1, 
    dif_rou ~ cab_seat_share + reld1,
    data = df_sd
)
msd2 <- selection(
    cabinet_party ~ par_seat_share + country_name + absd1 - 1, 
    dif_rou ~ cab_seat_share + absd1,
    data = df_sd
)
msd3 <- selection(
    cabinet_party ~ par_seat_share + country_name + reld196 - 1, 
    dif_rou ~ cab_seat_share + reld196,
    data = df_sd
)
msd4 <- selection(
    cabinet_party ~ par_seat_share + country_name + absd196 - 1, 
    dif_rou ~ cab_seat_share + absd196,
    data = df_sd
)


modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + rl - 1,
                make_matrix(df_sd, 3, c("exp", "ps", "rl"))
            ),
            msd1,
            c(1,35),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + rl + rho,
                make_matrix(dfc, 4, c("exp", "cs", "rl", "rho"))
            ),
            msd1,
            c(36:38),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ab - 1,
                make_matrix(df_sd, 3, c("exp", "ps", "ab"))
            ),
            msd2,
            c(1,35),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ab + rho,
                make_matrix(dfc, 4, c("exp", "cs", "ab", "rho"))
            ),
            msd2,
            c(36:38),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + rl - 1,
                make_matrix(df_sd, 3, c("exp", "ps", "rl"))
            ),
            msd3,
            c(1,35),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + rl + rho,
                make_matrix(dfc, 4, c("exp", "cs", "rl", "rho"))
            ),
            msd3,
            c(36:38),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ab - 1,
                make_matrix(df_sd, 3, c("exp", "ps", "ab"))
            ),
            msd4,
            c(1,35),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ab + rho,
                make_matrix(dfc, 4, c("exp", "cs", "ab", "rho"))
            ),
            msd4,
            c(36:38),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "rl" = "Relative radicals",
        "ab" = "Absolute radicals",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(msd1, 15, 17),
        heck_gof(msd2, 15, 17) |> select(-rowname),
        heck_gof(msd3, 15, 17) |> select(-rowname),
        heck_gof(msd4, 15, 17) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
add_header_above(
    c(
        "", 
        "$x_{ik}\\\\notin(\\\\bar{x_{k}}-\\\\sigma_{x_{k}}, \\\\bar{x_{k}}+\\\\sigma_{x_{k}})$" = 4, 
        "$x_{i}\\\\notin(\\\\bar{x_{ik}}-1.96\\\\sigma_{x_{k}}, \\\\bar{x_{k}}+1.96\\\\sigma_{x_{k}})$" = 4
    ), 
    escape = FALSE, 
    bold = TRUE
) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(10, hline_after = hl, extra_css = css) |>
row_spec(c(13:16), align = c("l", rep("c", 8)))
```

```{r}
#| label: tbl-fam
#| tbl-cap: Heckman models with party families
#| results: "asis"

mfam <- selection(
    cabinet_party ~ par_seat_share + country_name + family_name - 1, 
    dif_rou ~ cab_seat_share + family_name,
    data = df_fam
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + cd + so + co + ec + li + rw + sd + si - 1,
                make_matrix(df_fam, 10, c("exp", "ps", "cd", "so", "co", "ec", "li", "rw", "sd", "si"))
            ),
            mfam,
            c(1,35:42),
            type = "sel",
            cluster = df_fam$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + cd + so + co + ec + li + rw + sd + si + rho,
                make_matrix(df_fam_c, 11, c("exp", "cs", "cd", "so", "co", "ec", "li", "rw", "sd", "si", "rho"))
            ),
            mfam,
            c(43:52),
            type = "out",
            cluster = df_fam$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "so" = "Communist/Socialist",
        "ec" = "Green/Ecologist",
        "sd" = "Social democracy",
        "si" = "Special issue",
        "li" = "Liberal",
        "cd" = "Christian democracy",
        "co" = "Conservative",
        "rw" = "Right-wing",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mfam, 98, 99)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldd"
    } else if(knitr::is_html_output()){
        "lcc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
pack_rows("Party Family (Ref.: Agrarian)", 5, 20, italic = TRUE, bold = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(22, hline_after = hl, extra_css = css) |>
row_spec(c(25:28), align = c("l", rep("c", 2)))
```

```{r}
#| label: bargaining-data

dfbrg <- df |>
    filter(
        !is.na(gov_counter),
        !is.na(bargexp),
        !is.na(bargexp10),
        !is.na(bargexp60)
    )

dfbrg_c <- dfbrg |>
    filter(cabinet_party == 1)

dffp <- df |>
    filter(
        !is.na(prime_minister),
        !is.na(never_gov)
    )

dffp_c <- dffp |>
    filter(cabinet_party == 1)
```

```{r}
#| label: tbl-bargexp
#| tbl-cap: Heckman models with bargaining experience
#| results: "asis"

# TODO adjust percent coef labels

mbrg1 <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + gov_counter + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + gov_counter,
    data = dfbrg
)
mbrg2 <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + bargexp + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + bargexp,
    data = dfbrg
)
mbrg3 <-selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + bargexp10 + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + bargexp10,
    data = dfbrg
)
mbrg4 <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + bargexp60 + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + bargexp60,
    data = dfbrg
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + gc - 1,
                make_matrix(dfbrg, 4, c("exp", "ps", "pr", "gc"))
            ),
            mbrg1,
            c(1:3),
            type = "sel",
            cluster = dfbrg$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + gc + rho,
                make_matrix(dfbrg_c, 5, c("exp", "cs", "pr", "gc", "rho"))
            ),
            mbrg1,
            c(37:40),
            type = "out",
            cluster = dfbrg$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + bg - 1,
                make_matrix(dfbrg, 4, c("exp", "ps", "pr", "bg"))
            ),
            mbrg2,
            c(1:3),
            type = "sel",
            cluster = dfbrg$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + bg + rho,
                make_matrix(dfbrg_c, 5, c("exp", "cs", "pr", "bg", "rho"))
            ),
            mbrg2,
            c(37:40),
            type = "out",
            cluster = dfbrg$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + bg10 - 1,
                make_matrix(dfbrg, 4, c("exp", "ps", "pr", "bg10"))
            ),
            mbrg3,
            c(1:3),
            type = "sel",
            cluster = dfbrg$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + bg10 + rho,
                make_matrix(dfbrg_c, 5, c("exp", "cs", "pr", "bg10", "rho"))
            ),
            mbrg3,
            c(37:40),
            type = "out",
            cluster = dfbrg$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + bg60 - 1,
                make_matrix(dfbrg, 4, c("exp", "ps", "pr", "bg60"))
            ),
            mbrg4,
            c(1:3),
            type = "sel",
            cluster = dfbrg$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + bg60 + rho,
                make_matrix(dfbrg_c, 5, c("exp", "cs", "pr", "bg60", "rho"))
            ),
            mbrg4,
            c(37:40),
            type = "out",
            cluster = dfbrg$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr" = "Relative radicalism",
        "gc" = "Government counter",
        "bg" = "Bargaining experience (35$\\%$)",
        "bg10" = "Bargaining experience (10$\\%$)",
        "bg60" = "Bargaining experience (60$\\%$)",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mbrg1, 19, 21),
        heck_gof(mbrg2, 19, 21) |> select(-rowname),
        heck_gof(mbrg3, 19, 21) |> select(-rowname),
        heck_gof(mbrg4, 19, 21) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddddddd"
    } else if(knitr::is_html_output()){
        "lcccccccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2, "Model 4" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 8)))
```

```{r}
#| label: tbl-first-prime
#| tbl-cap: Heckman models with other controls
#| results: "asis"

mng <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + never_gov + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + never_gov,
    data = dffp
)
mpm <-selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + prime_minister + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + prime_minister,
    data = dffp
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + ng - 1,
                make_matrix(dffp, 4, c("exp", "ps", "pr", "ng"))
            ),
            mng,
            c(1:3),
            type = "sel",
            cluster = dffp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + ng + rho,
                make_matrix(dffp_c, 5, c("exp", "cs", "pr", "ng", "rho"))
            ),
            mng,
            c(37:40),
            type = "out",
            cluster = dffp$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + pm - 1,
                make_matrix(dffp, 4, c("exp", "ps", "pr", "pm"))
            ),
            mpm,
            c(1:3),
            type = "sel",
            cluster = dffp$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + pm + rho,
                make_matrix(dffp_c, 5, c("exp", "cs", "pr", "pm", "rho"))
            ),
            mpm,
            c(37:40),
            type = "out",
            cluster = dffp$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr" = "Relative radicalism",
        "ng" = "No previous government experience",
        "pm" = "Prime minister party",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mng, 15, 17),
        heck_gof(mpm, 15, 17) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddd"
    } else if(knitr::is_html_output()){
        "lcccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(12, hline_after = hl, extra_css = css) |>
row_spec(c(15:18), align = c("l", rep("c", 4)))
```

```{r}
#| label: tbl-no-seats
#| tbl-cap: Heckman models without seat share
#| results: "asis"

mst1 <- selection(
    cabinet_party ~ parlmean_dif_r + country_name - 1, 
    dif_rou ~ parlmean_dif_r,
    data = df
)
mst2 <- selection(
    cabinet_party ~ left_right_r + country_name - 1, 
    dif_rou ~ left_right_r,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ pr - 1,
                make_matrix(df, 2, c("exp", "pr"))
            ),
            mst1,
            1,
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ pr + rho,
                make_matrix(dfc, 3, c("exp", "pr", "rho"))
            ),
            mst1,
            c(35:36),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ lr - 1,
                make_matrix(df, 2, c("exp", "lr"))
            ),
            mst2,
            c(1),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ lr + rho,
                make_matrix(dfc, 3, c("exp", "lr", "rho"))
            ),
            mst2,
            c(35:36),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "pr" = "Relative radicalism",
        "lr" = "Absolute radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mst1, 9, 11),
        heck_gof(mst2, 9, 11) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        al <- "lcccc"
    } else if (knitr::is_latex_output()){
        al <- "ldddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(6, hline_after = hl, extra_css = css) |>
row_spec(c(9:12), align = c("l", rep("c", 4)))
```

```{r}
#| label: tbl-unrounded
#| tbl-cap: Heckman models with unrounded DV
#| results: "asis"

mur1 <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1, 
    dif_dec ~ cab_seat_share + parlmean_dif_r,
    data = df
)
mur2 <- selection(
    cabinet_party ~ par_seat_share + left_right_r + country_name - 1, 
    dif_dec ~ cab_seat_share + left_right_r,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr - 1,
                make_matrix(df, 3, c("exp", "ps", "pr"))
            ),
            mur1,
            c(1,2),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + rho,
                make_matrix(dfc, 4, c("exp", "cs", "pr", "rho"))
            ),
            mur1,
            c(36:38),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr - 1,
                make_matrix(df, 3, c("exp", "ps", "lr"))
            ),
            mur2,
            c(1,2),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + rho,
                make_matrix(dfc, 4, c("exp", "cs", "lr", "rho"))
            ),
            mur2,
            c(36:38),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr" = "Relative radicalism",
        "lr" = "Absolute radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mur1, 13, 15),
        heck_gof(mur2, 13, 15) |>
        select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldddd"
    } else if(knitr::is_html_output()){
        "lcccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(10, hline_after = hl, extra_css = css) |>
row_spec(c(13:16), align = c("l", rep("c", 4)))
```

```{r}
#| label: tbl-auth-cw-r
#| tbl-cap: Heckman models with interactions of relative radicalism with authoritarian legacy and pre-/post-Cold War status
#| results: "asis"

mauthr <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):auth_history + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):auth_history,
    data = df
)
mcwr <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):cold_war + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):cold_war,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1l + pr2l + pr1r + pr2r + pr1n + pr2n - 1,
                make_matrix(df, 8, c("exp", "ps", "pr1l", "pr2l", "pr1r", "pr2r", "pr1n", "pr2n"))
            ),
            mauthr,
            c(1,35:40),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1l + pr2l + pr1r + pr2r + pr1n + pr2n + rho,
                make_matrix(dfc, 9, c("exp", "cs", "pr1l", "pr2l", "pr1r", "pr2r", "pr1n", "pr2n", "rho"))
            ),
            mauthr,
            c(41:48),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1pr + pr2pr + pr1po + pr2po - 1,
                make_matrix(df, 6, c("exp", "ps", "pr1pr", "pr2pr", "pr1po", "pr2po"))
            ),
            mcwr,
            c(1,35:38),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1pr + pr2pr + pr1po + pr2po + rho,
                make_matrix(dfc, 7, c("exp", "cs", "pr1pr", "pr2pr", "pr1po", "pr2po", "rho"))
            ),
            mcwr,
            c(39:44),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr1l" = "Relative radicalism$\\times$Left Auth.",
        "pr2l" = "Relative radicalism$^{2}\\times$Left Auth.",
        "pr1r" = "Relative radicalism$\\times$Right Auth.",
        "pr2r" = "Relative radicalism$^{2}\\times$Right Auth.",
        "pr1n" = "Relative radicalism$\\times$Neither",
        "pr2n" = "Relative radicalism$^{2}\\times$Neither",
        "pr1pr" = "Relative radicalism$\\times$Pre Cold War end",
        "pr2pr" = "Relative radicalism$^{2}\\times$Pre Cold War end",
        "pr1po" = "Relative radicalism$\\times$Post Cold War end",
        "pr2po" = "Relative radicalism$^{2}\\times$Post Cold War end",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mauthr, 31, 33),
        heck_gof(mcwr, 31, 33) |>
        select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcccc"
    } else if (knitr::is_latex_output()){
        "ldddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(26, hline_after = hl, extra_css = css) |>
row_spec(c(29:32), align = c("l", rep("c", 4)))
```

```{r}
#| label: fig-auth-cw-r
#| fig-cap: "Predicted ministry compensation for relative radicalism interacted with authoritarian legacy and pre-/post-Cold War status. **Note**: Shaded area shows 95% confidence interval of this estimation. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1

ggarrange(
    predictions(
        mauthr,
        newdata = datagrid(
            parlmean_dif = rep(seq(-5, 5, .1), 3),
            auth_history = c(
                rep("right", 101),
                rep("left", 101),
                rep("neither", 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(mauthr, cluster = df$country_name)
    ) |>
    rename(int = auth_history) |>
    mutate(type = "Authoritarian History") |>
    ggplot(aes(parlmean_dif, estimate, colour = int, fill = int)) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = parlmean_dif, colour = auth_history), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "Predicted ministry compensation") +
        scale_x_continuous(
            name = "Relative Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        ),
    predictions(
        mcwr,
        newdata = datagrid(
            parlmean_dif = rep(seq(-5, 5, .1), 3),
            cold_war = c(
                rep("pre", 101),
                rep("post", 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(mcwr, cluster = df$country_name)
    ) |>
    rename(int = cold_war) |>
    mutate(type = "End of Cold War") |>
    ggplot(aes(parlmean_dif, estimate, colour = int, fill = int)) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = parlmean_dif, colour = cold_war), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "") +
        scale_x_continuous(
            name = "Relative Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        )
)
```

```{r}
#| label: tbl-auth-cw-a
#| tbl-cap: Heckman models with interactions of relative radicalism with authoritarian legacy and pre-/post-Cold War status
#| results: "asis"

mautha <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):auth_history + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):auth_history,
    data = df
)
mcwa <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):cold_war + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):cold_war,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1l + pr2l + pr1r + pr2r + pr1n + pr2n - 1,
                make_matrix(df, 8, c("exp", "ps", "pr1l", "pr2l", "pr1r", "pr2r", "pr1n", "pr2n"))
            ),
            mautha,
            c(1,35:40),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1l + pr2l + pr1r + pr2r + pr1n + pr2n + rho,
                make_matrix(dfc, 9, c("exp", "cs", "pr1l", "pr2l", "pr1r", "pr2r", "pr1n", "pr2n", "rho"))
            ),
            mautha,
            c(41:48),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1pr + pr2pr + pr1po + pr2po - 1,
                make_matrix(df, 6, c("exp", "ps", "pr1pr", "pr2pr", "pr1po", "pr2po"))
            ),
            mcwa,
            c(1,35:38),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1pr + pr2pr + pr1po + pr2po + rho,
                make_matrix(dfc, 7, c("exp", "cs", "pr1pr", "pr2pr", "pr1po", "pr2po", "rho"))
            ),
            mcwa,
            c(39:44),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr1l" = "Relative radicalism$\\times$Left Auth.",
        "pr2l" = "Relative radicalism$^{2}\\times$Left Auth.",
        "pr1r" = "Relative radicalism$\\times$Right Auth.",
        "pr2r" = "Relative radicalism$^{2}\\times$Right Auth.",
        "pr1n" = "Relative radicalism$\\times$Neither",
        "pr2n" = "Relative radicalism$^{2}\\times$Neither",
        "pr1pr" = "Relative radicalism$\\times$Pre Cold War end",
        "pr2pr" = "Relative radicalism$^{2}\\times$Pre Cold War end",
        "pr1po" = "Relative radicalism$\\times$Post Cold War end",
        "pr2po" = "Relative radicalism$^{2}\\times$Post Cold War end",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mautha, 31, 33),
        heck_gof(mcwa, 31, 33) |>
        select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcccc"
    } else if (knitr::is_latex_output()){
        "ldddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(26, hline_after = hl, extra_css = css) |>
row_spec(c(29:32), align = c("l", rep("c", 4)))
```

```{r}
#| label: fig-auth-cw-a
#| fig-cap: "Predicted ministry compensation for absolute radicalism interacted with authoritarian legacy and pre-/post-Cold War status. **Note**: Shaded area shows 95% confidence interval of this estimation. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1

ggarrange(
    predictions(
        mautha,
        newdata = datagrid(
            left_right = rep(seq(0, 10, .1), 3),
            auth_history = c(
                rep("right", 101),
                rep("left", 101),
                rep("neither", 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(mautha, cluster = df$country_name)
    ) |>
    rename(int = auth_history) |>
    mutate(type = "Authoritarian History") |>
    ggplot(aes(left_right, estimate, colour = int, fill = int)) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = left_right, colour = auth_history), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "Predicted ministry compensation") +
        scale_x_continuous(
            name = "Absolute Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        ),
    predictions(
        mcwa,
        newdata = datagrid(
            left_right = rep(seq(0, 10, .1), 3),
            cold_war = c(
                rep("pre", 101),
                rep("post", 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(mcwa, cluster = df$country_name)
    ) |>
    rename(int = cold_war) |>
    mutate(type = "End of Cold War") |>
    ggplot(aes(left_right, estimate, colour = int, fill = int)) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = left_right, colour = cold_war), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "") +
        scale_x_continuous(
            name = "Absolute Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        )
)
```

```{r}
#| label: tbl-time-r
#| tbl-cap: Heckman models with interactions of relative radicalism with time effects
#| results: "asis"

mtr1 <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):year + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):year,
    data = df
)
mtr2 <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):time_per + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):time_per,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1y + pr2y - 1,
                make_matrix(df, 4, c("exp", "ps", "pr1y", "pr2y"))
            ),
            mtr1,
            c(1,35, 36),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1y + pr2y + rho,
                make_matrix(dfc, 5, c("exp", "cs", "pr1y", "pr2y", "rho"))
            ),
            mtr1,
            c(37:40),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1a + pr2a + pr1b + pr2b + pr1c + pr2c + pr1d + pr2d + pr1e + pr2e - 1,
                make_matrix(df, 12, c("exp", "ps", "pr1a", "pr2a", "pr1b", "pr2b", "pr1c", "pr2c", "pr1d", "pr2d", "pr1e", "pr2e"))
            ),
            mtr2,
            c(1,35:44),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1a + pr2a + pr1b + pr2b + pr1c + pr2c + pr1d + pr2d + pr1e + pr2e + rho,
                make_matrix(dfc,13, c("exp", "cs", "pr1a", "pr2a", "pr1b", "pr2b", "pr1c", "pr2c", "pr1d", "pr2d", "pr1e", "pr2e", "rho"))
            ),
            mtr2,
            c(45:56),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr1y" = "Relative radicalism$\\times$Year",
        "pr2y" = "Relative radicalism$^{2}\\times$Year",
        "pr1a" = "Relative radicalism$\\times$ 1960 and earlier",
        "pr2a" = "Relative radicalism$^{2}\\times$ 1960 and earlier",
        "pr1b" = "Relative radicalism$\\times$ 1961 to 1980",
        "pr2b" = "Relative radicalism$^{2}\\times$ 1961 to 1980",
        "pr1c" = "Relative radicalism$\\times$ 1981 to 2000",
        "pr2c" = "Relative radicalism$^{2}\\times$ 1981 to 2000",
        "pr1d" = "Relative radicalism$\\times$ 2001 and 2020",
        "pr2d" = "Relative radicalism$^{2}\\times$ 2001 and 2020",
        "pr1e" = "Relative radicalism$\\times$ 2021 and later",
        "pr2e" = "Relative radicalism$^{2}\\times$ 2021 and later",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mtr1, 35, 37),
        heck_gof(mtr2, 35, 37) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcccc"
    } else if (knitr::is_latex_output()){
        "ldddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(30, hline_after = hl, extra_css = css) |>
row_spec(c(33:36), align = c("l", rep("c", 4)))
```

```{r}
#| label: fig-time-r
#| fig-cap: "Predicted ministry compensation for relative radicalism interacted with time effects. **Note**: Shaded area shows 95% confidence interval of this estimation. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1

ggarrange(
    predictions(
        mtr1,
        newdata = datagrid(
            parlmean_dif = rep(seq(-5, 5, .1), 3),
            year = c(
                rep(1960, 101),
                rep(1990, 101),
                rep(2020, 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(mtr1, cluster = df$country_name)
    ) |>
    rename(int = year) |>
    mutate(type = "Linear Time Trend") |>
    ggplot(aes(parlmean_dif, estimate, colour = as.factor(int), fill = as.factor(int))) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = parlmean_dif), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "Predicted ministry compensation") +
        scale_x_continuous(
            name = "Relative Radicalism",
            expand = c(0,0)
        ) +
        guides(
            fill=guide_legend(ncol=2),
            colour=guide_legend(ncol=2)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        ),
    predictions(
        mtr2,
        newdata = datagrid(
            parlmean_dif = rep(seq(-5, 5, .1), 4),
            time_per = c(
                rep("1960 and earlier", 101),
                rep("1961-1980", 101),
                rep("1981-2000", 101),
                rep("2001-2020", 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(mtr2, cluster = df$country_name)
    ) |>
    rename(int = time_per) |>
    mutate(type = "Time Period") |>
    ggplot(aes(parlmean_dif, estimate, colour = int, fill = int)) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc |> filter(time_per != "2021 and later"), aes(x = parlmean_dif, colour = time_per), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "") +
        scale_x_continuous(
            name = "Relative Radicalism",
            expand = c(0,0)
        ) +
        guides(
            fill=guide_legend(ncol=2),
            colour=guide_legend(ncol=2)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        )
)
```

```{r}
#| label: tbl-time-a
#| tbl-cap: Heckman models with interactions of absolute radicalism with time period effect
#| results: "asis"

mta <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):time_per + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):time_per,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1a + pr2a + pr1b + pr2b + pr1c + pr2c + pr1d + pr2d + pr1e + pr2e - 1,
                make_matrix(df, 12, c("exp", "ps", "pr1a", "pr2a", "pr1b", "pr2b", "pr1c", "pr2c", "pr1d", "pr2d", "pr1e", "pr2e"))
            ),
            mta,
            c(1,35:44),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1a + pr2a + pr1b + pr2b + pr1c + pr2c + pr1d + pr2d + pr1e + pr2e + rho,
                make_matrix(dfc,13, c("exp", "cs", "pr1a", "pr2a", "pr1b", "pr2b", "pr1c", "pr2c", "pr1d", "pr2d", "pr1e", "pr2e", "rho"))
            ),
            mta,
            c(45:56),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr1a" = "Absolute radicalism$\\times$ 1960 and earlier",
        "pr2a" = "Absolute radicalism$^{2}\\times$ 1960 and earlier",
        "pr1b" = "Absolute radicalism$\\times$ 1961 to 1980",
        "pr2b" = "Absolute radicalism$^{2}\\times$ 1961 to 1980",
        "pr1c" = "Absolute radicalism$\\times$ 1981 to 2000",
        "pr2c" = "Absolute radicalism$^{2}\\times$ 1981 to 2000",
        "pr1d" = "Absolute radicalism$\\times$ 2001 and 2020",
        "pr2d" = "Absolute radicalism$^{2}\\times$ 2001 and 2020",
        "pr1e" = "Absolute radicalism$\\times$ 2021 and later",
        "pr2e" = "Absolute radicalism$^{2}\\times$ 2021 and later",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(heck_gof(mta, 31, 33)),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcc"
    } else if (knitr::is_latex_output()){
        "ldd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(26, hline_after = hl, extra_css = css) |>
row_spec(c(29:32), align = c("l", rep("c", 2)))
```

```{r}
#| label: fig-time-a
#| fig-cap: "Predicted ministry compensation for absolute radicalism interacted with time period. **Note**: Shaded area shows 95% confidence interval of this estimation. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1

predictions(
    mta,
    newdata = datagrid(
        left_right = rep(seq(0, 10, .1), 3),
        time_per = c(
            rep("1960 and earlier", 101),
            rep("1961-1980", 101),
            rep("1981-2000", 101),
            rep("2001-2020", 101)
        )
    ),
    type = "unconditional",
    vcov = vcovCL(mta, cluster = df$country_name)
) |>
rename(int = time_per) |>
mutate(type = "Time Period") |>
ggplot(aes(left_right, estimate, colour = int, fill = int)) +
    facet_wrap(~type) +
    geom_line() +
    geom_rug(data = dfc |> filter(time_per != "2021 and later"), aes(x = left_right, colour = time_per), inherit.aes = FALSE, alpha = .2) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
    scale_y_continuous(name = "Predicted ministry compensation") +
    scale_x_continuous(
        name = "Absolute Radicalism",
        expand = c(0,0)
    ) +
    theme(
        legend.title = element_blank(),
        legend.position = "bottom"
    )
```

```{r}
#| label: tbl-eu-r
#| tbl-cap: Heckman models with interactions of relative radicalism with EU effects
#| results: "asis"

meur1 <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):eu_member + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):eu_member,
    data = df
)
meur2 <- selection(
    cabinet_party ~ par_seat_share + poly(parlmean_dif, 2, raw = TRUE):seats_eu + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(parlmean_dif, 2, raw = TRUE):seats_eu,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1e + pr2e + pr1n + pr2n - 1,
                make_matrix(df, 6, c("exp", "ps", "pr1e", "pr2e", "pr1n", "pr2n"))
            ),
            meur1,
            c(1,35:38),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1e + pr2e + pr1n + pr2n + rho,
                make_matrix(dfc, 7, c("exp", "cs", "pr1e", "pr2e", "pr1n", "pr2n", "rho"))
            ),
            meur1,
            c(39:44),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1r + pr2r + pr1x + pr2x - 1,
                make_matrix(df, 6, c("exp", "ps", "pr1r", "pr2r", "pr1x", "pr2x"))
            ),
            meur2,
            c(1,35:38),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1r + pr2r + pr1x + pr2x + rho,
                make_matrix(dfc, 7, c("exp", "cs", "pr1r", "pr2r", "pr1x", "pr2x", "rho"))
            ),
            meur2,
            c(39:44),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr1e" = "Relative radicalism$\\times$EU member state",
        "pr2e" = "Relative radicalism$^{2}\\times$EU member state",
        "pr1n" = "Relative radicalism$\\times$Non-EU state",
        "pr2n" = "Relative radicalism$^{2}\\times$Non-EU state",
        "pr1r" = "Relative radicalism$\\times$In EP",
        "pr2r" = "Relative radicalism$^{2}\\times$In EP",
        "pr1x" = "Relative radicalism$\\times$Not in EP",
        "pr2x" = "Relative radicalism$^{2}\\times$Not in EP",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(meur1, 27, 29),
        heck_gof(meur2, 27, 29) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcccc"
    } else if (knitr::is_latex_output()){
        "ldddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(22, hline_after = hl, extra_css = css) |>
row_spec(c(25:28), align = c("l", rep("c", 4)))
```

```{r}
#| label: fig-eu-r
#| fig-cap: "Predicted ministry compensation for relative radicalism interacted with EU effects. **Note**: Shaded area shows 95% confidence interval of this estimation. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1

ggarrange(
    predictions(
        meur1,
        newdata = datagrid(
            parlmean_dif = rep(seq(-5, 5, .1), 2),
            eu_member = c(
                rep(TRUE, 101),
                rep(FALSE, 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(meur1, cluster = df$country_name)
    ) |>
    rename(int = eu_member) |>
    mutate(type = "EU Membership") |>
    ggplot(aes(parlmean_dif, estimate, colour = as.factor(int), fill = as.factor(int))) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = parlmean_dif, colour = eu_member), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "Predicted ministry compensation") +
        scale_x_continuous(
            name = "Relative Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        ),
    predictions(
        meur2,
        newdata = datagrid(
            parlmean_dif = rep(seq(-5, 5, .1), 2),
            seats_eu = c(
                rep(TRUE, 101),
                rep(FALSE, 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(meur2, cluster = df$country_name)
    ) |>
    rename(int = seats_eu) |>
    mutate(type = "Party in European Parliament") |>
    ggplot(aes(parlmean_dif, estimate, colour = int, fill = int)) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = parlmean_dif, colour = seats_eu), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "") +
        scale_x_continuous(
            name = "Relative Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        )
)
```

```{r}
#| label: tbl-eu-a
#| tbl-cap: Heckman models with interactions of absolute radicalism with EU effects
#| results: "asis"

meua1 <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):eu_member + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):eu_member,
    data = df
)
meua2 <- selection(
    cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):seats_eu + country_name - 1, 
    dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):seats_eu,
    data = df
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1e + pr2e + pr1n + pr2n - 1,
                make_matrix(df, 6, c("exp", "ps", "pr1e", "pr2e", "pr1n", "pr2n"))
            ),
            meua1,
            c(1,35:38),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1e + pr2e + pr1n + pr2n + rho,
                make_matrix(dfc, 7, c("exp", "cs", "pr1e", "pr2e", "pr1n", "pr2n", "rho"))
            ),
            meua1,
            c(39:44),
            type = "out"
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr1r + pr2r + pr1x + pr2x - 1,
                make_matrix(df, 6, c("exp", "ps", "pr1r", "pr2r", "pr1x", "pr2x"))
            ),
            meua2,
            c(1,35:38),
            type = "sel"
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr1r + pr2r + pr1x + pr2x + rho,
                make_matrix(dfc, 7, c("exp", "cs", "pr1r", "pr2r", "pr1x", "pr2x", "rho"))
            ),
            meua2,
            c(39:44),
            type = "out"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr1e" = "Absolute radicalism$\\times$EU member state",
        "pr2e" = "Absolute radicalism$^{2}\\times$EU member state",
        "pr1n" = "Absolute radicalism$\\times$Non-EU state",
        "pr2n" = "Absolute radicalism$^{2}\\times$Non-EU state",
        "pr1r" = "Absolute radicalism$\\times$In EP",
        "pr2r" = "Absolute radicalism$^{2}\\times$In EP",
        "pr1x" = "Absolute radicalism$\\times$Not in EP",
        "pr2x" = "Absolute radicalism$^{2}\\times$Not in EP",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(meua1, 27, 29),
        heck_gof(meua2, 27, 29) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcccc"
    } else if (knitr::is_latex_output()){
        "ldddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(22, hline_after = hl, extra_css = css) |>
row_spec(c(25:28), align = c("l", rep("c", 4)))
```

```{r}
#| label: fig-eu-a
#| fig-cap: "Predicted ministry compensation for absolute radicalism interacted with EU effects. **Note**: Shaded area shows 95% confidence interval of this estimation. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1

ggarrange(
    predictions(
        meua1,
        newdata = datagrid(
            left_right = rep(seq(0, 10, .1), 2),
            eu_member = c(
                rep(TRUE, 101),
                rep(FALSE, 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(meua1, cluster = df$country_name)
    ) |>
    rename(int = eu_member) |>
    mutate(type = "EU Membership") |>
    ggplot(aes(left_right, estimate, colour = as.factor(int), fill = as.factor(int))) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = left_right, colour = eu_member), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "Predicted ministry compensation") +
        scale_x_continuous(
            name = "Absolute Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        ),
    predictions(
        meua2,
        newdata = datagrid(
            left_right = rep(seq(0, 10, .1), 2),
            seats_eu = c(
                rep(TRUE, 101),
                rep(FALSE, 101)
            )
        ),
        type = "unconditional",
        vcov = vcovCL(meua1, cluster = df$country_name)
    ) |>
    rename(int = seats_eu) |>
    mutate(type = "Party in European Parliament") |>
    ggplot(aes(left_right, estimate, colour = int, fill = int)) +
        facet_wrap(~type) +
        geom_line() +
        geom_rug(data = dfc, aes(x = left_right, colour = seats_eu), inherit.aes = FALSE, alpha = .2) +
        geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, colour = NA) +
        scale_y_continuous(name = "") +
        scale_x_continuous(
            name = "Absolute Radicalism",
            expand = c(0,0)
        ) +
        theme(
            legend.title = element_blank(),
            legend.position = "bottom"
        )
)
```

```{r}
#| label: interaction-prep

pr_cn <- predictions(
        selection(
            cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):country_name,
            dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):country_name,
            data = df
        ),
        data.frame(cn = sort(rep(unique(na.omit(df$country_name)), 101))) |>
            mutate(
                country_name = cn,
                left_right = rep(seq(0,10,.1), 33),
                cab_seat_share = .25
            ) |>
            select(-cn),
        type = "unconditional",
        vcov = vcovCL(
            selection(
                cabinet_party ~ par_seat_share + poly(left_right, 2, raw = TRUE):country_name,
                dif_rou ~ cab_seat_share + poly(left_right, 2, raw = TRUE):country_name,
                data = df
            ),
            cluster = df$country_name
        )
    ) |>
    left_join(
        df |>
            select(country_name, country_name_short),
        by = "country_name"
    ) |>
    distinct() |>
    drop_na(country_name_short)

pr_cn_max <- pr_cn |>
    group_by(country_name_short) |>
    mutate(lr_max = left_right[estimate == max(estimate)]) |>
    summarise(lr_max = first(lr_max))

nou <- sum(pr_cn_max$lr_max == 10 | pr_cn_max$lr_max == 0)
cou <- n_distinct(pr_cn_max$country_name_short)
```

```{r}
#| label: fig-interaction
#| fig-cap: 'Country Interactions. **Note**: Predicted deviation from a proportional allocation of ministries per country. 95 per cent confidence interval displayed as ribbon.'
#| fig-height: 7
#| fig-width: 6.1

ggplot(pr_cn, aes(left_right, estimate)) +
    facet_wrap(~country_name_short, ncol = 5, scale = "free_y") +
    geom_hline(yintercept = 0, colour = "grey20", linetype = "dashed") +
    geom_vline(xintercept = 5, colour = "grey20", linetype = "dashed") +
    geom_vline(
        data = pr_cn_max,
        aes(xintercept = lr_max, colour = lr_max),
        linetype = "solid"
    ) +
    geom_line() +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3) +
    scale_colour_gradient(
        name = "",
        limits = c(0,10),
        high = "mediumblue",
        low = "red3",
        guide = "none"
    ) +
    scale_y_continuous(name = "Deviation") +
    scale_x_continuous(
        name = "Left-right position",
        breaks = seq(0,10,2.5),
        expand = c(0,0),
        labels = c("left", "", "centre", "", "right")
    ) +
    theme(
        panel.grid.minor = element_blank(),
        plot.margin = margin(2, 4, 2, 2, unit = "mm")
    )
```

```{r}
#| label: tbl-countries1
#| tbl-cap: Heckman models excluding AUS, AUT, BEL, and BGR
#| results: "asis"

cn_exl_tbl(c("AUS", "AUT", "BEL", "BGR")) |>
    add_header_above(
        c(
            "",
            "No AUS" = 2,
            "No AUT" = 2,
            "No BEL" = 2,
            "No BGR" = 2
        ),
        bold = TRUE
    )
```

```{r}
#| label: tbl-countries2
#| tbl-cap: Heckman models excluding CHE, CZE, DEU, and DNK
#| results: "asis"

cn_exl_tbl(c("CHE", "CZE", "DEU", "DNK")) |>
    add_header_above(
        c(
            "",
            "No CHE" = 2,
            "No CZE" = 2,
            "No DEU" = 2,
            "No DNK" = 2
        ),
        bold = TRUE
    )
```

```{r}
#| label: tbl-countries3
#| tbl-cap: Heckman models excluding ESP, EST, FIN, and FRA
#| results: "asis"

cn_exl_tbl(c("ESP", "EST", "FIN", "FRA")) |>
    add_header_above(
        c(
            "",
            "No ESP" = 2,
            "No EST" = 2,
            "No FIN" = 2,
            "No FRA" = 2
        ),
        bold = TRUE
    )
```

```{r}
#| label: tbl-countries4
#| tbl-cap: Heckman models excluding GBR, GRC, HRV, and HUN
#| results: "asis"

cn_exl_tbl(c("GBR", "GRC", "HRV", "HUN")) |>
    add_header_above(
        c(
            "",
            "No GBR" = 2,
            "No GRC" = 2,
            "No HRV" = 2,
            "No HUN" = 2
        ),
        bold = TRUE
    )
```

```{r}
#| label: tbl-countries5
#| tbl-cap: Heckman models excluding IRL, ISL, ISR, and ITA
#| results: "asis"

cn_exl_tbl(c("IRL", "ISL", "ISR", "ITA")) |>
    add_header_above(
        c(
            "",
            "No IRL" = 2,
            "No ISL" = 2,
            "No ISR" = 2,
            "No ITA" = 2
        ),
        bold = TRUE
    )
```

```{r}
#| label: tbl-countries6
#| tbl-cap: Heckman models excluding JPN, LTU, LUX, and LVA
#| results: "asis"

cn_exl_tbl(c("JPN", "LTU", "LUX", "LVA")) |>
    add_header_above(
        c(
            "",
            "No JPN" = 2,
            "No LTU" = 2,
            "No LUX" = 2,
            "No LVA" = 2
        ),
        bold = TRUE
    )
```

```{r}
#| label: tbl-countries7
#| tbl-cap: Heckman models excluding NLD, NOR, NZL, and POL
#| results: "asis"

cn_exl_tbl(c("NLD", "NOR", "NZL", "POL")) |>
    add_header_above(
        c(
            "",
            "No NLD" = 2,
            "No NOR" = 2,
            "No NZL" = 2,
            "No POL" = 2
        ),
        bold = TRUE
    )
```

```{r}
#| label: tbl-countries8
#| tbl-cap: Heckman models excluding PRT, ROU, SVK, and SVN
#| results: "asis"

cn_exl_tbl(c("PRT", "ROU", "SVK", "SVN")) |>
    add_header_above(
        c(
            "",
            "No PRT" = 2,
            "No ROU" = 2,
            "No SVK" = 2,
            "No SVN" = 2
        ),
        bold = TRUE
    )
```

\elandscape

```{r}
#| label: tbl-countries9
#| tbl-cap: Heckman model excluding SWE
#| results: "asis"

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr - 1,
                make_matrix(
                    df |> filter(country_name_short != "SWE"),
                    3,
                    c("exp", "ps", "pr")
                )
            ),
            selection(
                cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                dif_rou ~ cab_seat_share + parlmean_dif_r,
                data = df |> filter(country_name_short != "SWE")
            ),
            c(1:2),
            type = "sel",
            cluster = df |> 
                filter(country_name_short != "SWE") |> 
                select(country_name)
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + rho,
                make_matrix(
                    df |> filter(country_name_short != "SWE"),
                    4,
                    c("exp", "cs", "pr", "rho")
                )
            ),
            selection(
                cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
                dif_rou ~ cab_seat_share + parlmean_dif_r,
                data = df |> filter(country_name_short != "SWE")
            ),
            c(35:37),
            type = "out",
            cluster = df |> 
                filter(country_name_short != "SWE") |> 
                select(country_name)
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr" = "Relative radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = heck_gof(
        selection(
            cabinet_party ~ par_seat_share + parlmean_dif_r + country_name - 1,
            dif_rou ~ cab_seat_share + parlmean_dif_r,
            data = df |> filter(country_name_short != "SWE")
        ), 12, 13
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if(knitr::is_latex_output()){
        "ldd"
    } else if(knitr::is_html_output()){
        "lcc"
    },
    output = "latex",
    escape = FALSE
) |>
kable_styling(latex_options = ltop, full_width = TRUE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(8, hline_after = hl, extra_css = css) |>
row_spec(c(11:14), align = c("l", rep("c", 2))) |>
add_header_above(c("", "No SWE" = 2), bold = TRUE)
```

```{r}
#| label: tbl-sel-dif-r
#| tbl-cap: Heckman models accounting for relative radicalisation
#| results: "asis"

df_dif <- df |>
    filter(!is.na(parlmean_dif_r_dif))

df_dif_c <- df_dif |>
    filter(cabinet_party == 1)

mcsr <- selection(
    cabinet_party ~ par_seat_share + parlmean_dif_r + parlmean_dif_r_dif + country_name - 1, 
    dif_rou ~ cab_seat_share + parlmean_dif_r + parlmean_dif_r_dif,
    data = df_dif
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + pr + prd - 1,
                make_matrix(df_dif, 4, c("exp", "ps", "pr", "prd"))
            ),
            mcsr,
            c(1:3),
            type = "sel",
            cluster = df_dif$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + pr + prd + rho,
                make_matrix(df_dif_c, 5, c("exp", "cs", "pr", "prd", "rho"))
            ),
            mcsr,
            c(37:40),
            type = "out",
            cluster = df_dif$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "pr" = "Relative radicalism",
        "prd" = "$\\Delta$Relative Radicalism",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(heck_gof(mcsr, 15, 17)),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcc"
    } else if (knitr::is_latex_output()){
        "ldd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(10, hline_after = hl, extra_css = css) |>
row_spec(c(13:16), align = c("l", rep("c", 2)))
```

```{r}
#| label: tbl-sel-dif-ches
#| tbl-cap: Heckman models accounting for CHES-based absolute radicalisation
#| results: "asis"

df_dif <- df |>
    filter(
        !is.na(lrgen_r_dif), 
        !is.na(lrecon_r_dif),
        !is.na(galtan_r_dif)
    )

df_dif_c <- df_dif |>
    filter(cabinet_party == 1)

mcsches1 <- selection(
    cabinet_party ~ par_seat_share + lrgen_r + lrgen_r_dif + country_name - 1, 
    dif_rou ~ cab_seat_share + lrgen_r + lrgen_r_dif,
    data = df_dif
)
mcsches2 <- selection(
    cabinet_party ~ par_seat_share + lrecon_r + lrecon_r_dif + country_name - 1, 
    dif_rou ~ cab_seat_share + lrecon_r + lrecon_r_dif,
    data = df_dif
)
mcsches3 <- selection(
    cabinet_party ~ par_seat_share + galtan_r + galtan_r_dif + country_name - 1, 
    dif_rou ~ cab_seat_share + galtan_r + galtan_r_dif,
    data = df_dif
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr + lrd - 1,
                make_matrix(df_dif, 4, c("exp", "ps", "lr", "lrd"))
            ),
            mcsches1,
            c(1:3),
            type = "sel",
            cluster = df_dif$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + lrd + rho,
                make_matrix(df_dif_c, 5, c("exp", "cs", "lr", "lrd", "rho"))
            ),
            mcsches1,
            c(29:32),
            type = "out",
            cluster = df_dif$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + ecd - 1,
                make_matrix(df_dif, 4, c("exp", "ps", "ec", "ecd"))
            ),
            mcsches2,
            c(1:3),
            type = "sel",
            cluster = df_dif$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + ecd + rho,
                make_matrix(df_dif_c, 5, c("exp", "cs", "ec", "ecd", "rho"))
            ),
            mcsches2,
            c(29:32),
            type = "out",
            cluster = df_dif$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + gt + gtd - 1,
                make_matrix(df_dif, 4, c("exp", "ps", "gt", "gtd"))
            ),
            mcsches3,
            c(1:3),
            type = "sel",
            cluster = df_dif$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + gt + gtd + rho,
                make_matrix(df_dif_c, 5, c("exp", "cs", "gt", "gtd", "rho"))
            ),
            mcsches3,
            c(29:32),
            type = "out",
            cluster = df_dif$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "lr" = "Absolute radicalism (Left-Right)",
        "lrd" = "$\\Delta$Absolute radicalism (Left-Right)",
        "ec" = "Absolute radicalism (Economic)",
        "ecd" = "$\\Delta$Absolute radicalism (Economic)",
        "gt" = "Absolute radicalism (Cultural)",
        "gtd" = "$\\Delta$Absolute radicalism (Cultural)",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mcsches1, 23, 25),
        heck_gof(mcsches2, 23, 25) |> select(-rowname),
        heck_gof(mcsches3, 23, 25) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcccccc"
    } else if (knitr::is_latex_output()){
        "ldddddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(18, hline_after = hl, extra_css = css) |>
row_spec(c(21:24), align = c("l", rep("c", 6)))
```

```{r}
#| label: tbl-sel-dif-cmp
#| tbl-cap: Heckman models accounting for CMP-based absolute radicalisation
#| results: "asis"

df_dif <- df |>
    filter(
        !is.na(rile_r_dif), 
        !is.na(market_nat_r_dif),
        !is.na(culture_nat_r_dif)
    )

df_dif_c <- df_dif |>
    filter(cabinet_party == 1)

mcscmp1 <- selection(
    cabinet_party ~ par_seat_share + rile_r + rile_r_dif + country_name - 1, 
    dif_rou ~ cab_seat_share + rile_r + rile_r_dif,
    data = df_dif
)
mcscmp2 <- selection(
    cabinet_party ~ par_seat_share + market_nat_r + market_nat_r_dif + country_name - 1, 
    dif_rou ~ cab_seat_share + market_nat_r + market_nat_r_dif,
    data = df_dif
)
mcscmp3 <- selection(
    cabinet_party ~ par_seat_share + culture_nat_r + culture_nat_r_dif + country_name - 1, 
    dif_rou ~ cab_seat_share + culture_nat_r + culture_nat_r_dif,
    data = df_dif
)

modelsummary(
    list(
        "Selection" = change_coefs(
            lm(
                exp ~ ps + lr + lrd - 1,
                make_matrix(df_dif, 4, c("exp", "ps", "lr", "lrd"))
            ),
            mcscmp1,
            c(1:3),
            type = "sel",
            cluster = df_dif$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + lr + lrd + rho,
                make_matrix(df_dif_c, 5, c("exp", "cs", "lr", "lrd", "rho"))
            ),
            mcscmp1,
            c(37:40),
            type = "out",
            cluster = df_dif$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + ec + ecd - 1,
                make_matrix(df_dif, 4, c("exp", "ps", "ec", "ecd"))
            ),
            mcscmp2,
            c(1:3),
            type = "sel",
            cluster = df_dif$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + ec + ecd + rho,
                make_matrix(df_dif_c, 5, c("exp", "cs", "ec", "ecd", "rho"))
            ),
            mcscmp2,
            c(37:40),
            type = "out",
            cluster = df_dif$country_name
        ),
        "Selection" = change_coefs(
            lm(
                exp ~ ps + gt + gtd - 1,
                make_matrix(df_dif, 4, c("exp", "ps", "gt", "gtd"))
            ),
            mcscmp3,
            c(1:3),
            type = "sel",
            cluster = df_dif$country_name
        ),
        "Outcome" = change_coefs(
            lm(
                exp ~ cs + gt + gtd + rho,
                make_matrix(df_dif_c, 5, c("exp", "cs", "gt", "gtd", "rho"))
            ),
            mcscmp3,
            c(37:40),
            type = "out",
            cluster = df_dif$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "ps" = "Parliamentary seat share",
        "cs" = "Cabinet seat share",
        "lr" = "Absolute radicalism (Left-Right)",
        "lrd" = "$\\Delta$Absolute radicalism (Left-Right)",
        "ec" = "Absolute radicalism (Economic)",
        "ecd" = "$\\Delta$Absolute radicalism (Economic)",
        "gt" = "Absolute radicalism (Cultural)",
        "gtd" = "$\\Delta$Absolute radicalism (Cultural)",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(mcscmp1, 23, 25),
        heck_gof(mcscmp2, 23, 25) |> select(-rowname),
        heck_gof(mcscmp3, 23, 25) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcccccc"
    } else if (knitr::is_latex_output()){
        "ldddddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2, "Model 3" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(18, hline_after = hl, extra_css = css) |>
row_spec(c(21:24), align = c("l", rep("c", 6)))
```

```{r}
#| label: tbl-cau-cont
#| tbl-cap: Linear models for compensation change in continued coalitions
#| results: "asis"

dfcont <- df |>
    filter(
        continuation == TRUE,
        cabinet_party == 1
    )

modelsummary(
    list(
        "Model 1" = coeftest(
            lm(dif_rou_dif ~ cab_seat_share_dif + cab_seat_share_l + parlmean_dif_r_dif + parlmean_dif_rl + dif_rou_l, dfcont),
            vcovCL,
            cluster = dfcont$country_name
        ),
        "Model 2" = coeftest(
            lm(dif_rou_dif ~ cab_seat_share_dif + cab_seat_share_l + lrgen_r_dif + lrgen_rl + dif_rou_l, dfcont),
            vcovCL,
            cluster = dfcont$country_name
        ),
        "Model 3" = coeftest(
            lm(dif_rou_dif ~ cab_seat_share_dif + cab_seat_share_l + lrecon_r_dif + lrecon_rl + dif_rou_l, dfcont),
            vcovCL,
            cluster = dfcont$country_name
        ),
        "Model 4" = coeftest(
            lm(dif_rou_dif ~ cab_seat_share_dif + cab_seat_share_l + galtan_r_dif + galtan_rl + dif_rou_l, dfcont),
            vcovCL,
            cluster = dfcont$country_name
        ),
        "Model 5" = coeftest(
            lm(dif_rou_dif ~ cab_seat_share_dif + cab_seat_share_l + rile_r_dif + rile_rl + dif_rou_l, dfcont),
            vcovCL,
            cluster = dfcont$country_name
        ),
        "Model 6" = coeftest(
            lm(dif_rou_dif ~ cab_seat_share_dif + cab_seat_share_l + market_nat_r_dif + market_nat_rl + dif_rou_l, dfcont),
            vcovCL,
            cluster = dfcont$country_name
        ),
        "Model 7" = coeftest(
            lm(dif_rou_dif ~ cab_seat_share_dif + cab_seat_share_l + culture_nat_r_dif + culture_nat_rl + dif_rou_l, dfcont),
            vcovCL,
            cluster = dfcont$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "cab_seat_share_dif" = "Cabinet seat share",
        "parlmean_dif_rl" = "Relative Radicalism$_{t-1}$",
        "parlmean_dif_r_dif" = "$\\Delta$Relative Radicalism",
        "lrgen_rl" = "Absolute radicalism (Left-Right)$_{t-1}$",
        "rile_rl" = "Absolute radicalism (Left-Right)$_{t-1}$",
        "lrgen_r_dif" = "$\\Delta$Absolute radicalism (Left-Right)",
        "rile_r_dif" = "$\\Delta$Absolute radicalism (Left-Right)",
        "lrecon_rl" = "Absolute radicalism (Economic)$_{t-1}$",
        "market_nat_rl" = "Absolute radicalism (Economic)$_{t-1}$",
        "lrecon_r_dif" = "$\\Delta$Absolute radicalism (Economic)",
        "market_nat_r_dif" = "$\\Delta$Absolute radicalism (Economic)",
        "galtan_rl" = "Absolute radicalism (Cultural)$_{t-1}$",
        "culture_nat_rl" = "Absolute radicalism (Cultural)$_{t-1}$",
        "galtan_r_dif" = "$\\Delta$Absolute radicalism (Cultural)",
        "culture_nat_r_dif" = "$\\Delta$Absolute radicalism (Cultural)",
        "dif_rou_l" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept"
    ),
    gof_map = list(
        list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f"),
        list("raw" = "aic", "clean" = "AIC", "fmt" = "%.1f"),
        list("raw" = "bic", "clean" = "BIC", "fmt" = "%.1f")
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccccccc"
    } else if (knitr::is_latex_output()){
        "lddddddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "", "CHES" = 3, "Manifesto" = 3), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(c(23:25), align = c("l", rep("c", 7)))
```

```{r}
#| label: three-step-data

dfh <- df |>
    group_by(party_id) |>
    filter(
        !is.na(dif_rou_dif),
        previous_cabinet_id == lag(cabinet_id)
    ) |>
    ungroup()

dfhc <- dfh |> filter(cab_lag == 1)
```

```{r}
#| label: tbl-three-r
#| tbl-cap: Three-step Heckman models accounting for relative radicalisation
#| results: "asis"

sel <- selection(
        cab_lag ~ par_seat_share_l + parlmean_dif_rl + country_name - 1,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + parlmean_dif_rl + parlmean_dif_r_dif,
        data = dfh
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfh, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + parlmean_dif_rl + parlmean_dif_r_dif + dif_rou_l + imr,
    data = dfh2
)

modelsummary(
    list(
        "Initial Inclusion" = change_coefs(
            lm(
                exp ~ psl + prl - 1,
                make_matrix(dfh, 3, c("exp", "psl", "prl"))
            ),
            sel,
            c(1,2),
            type = "sel",
            cluster= dfh$country_name
        ),
        "Continued Inclusion" = change_coefs(
            lm(
                exp ~ psl + psd + prl + prd + rho,
                make_matrix(dfhc, 6, c("exp", "psl", "psd", "prl", "prd", "rho"))
            ),
            sel,
            c(34:38),
            type = "out",
            cluster= dfh$country_name
        ),
        "$\\Delta$ Gamson Deviation" = change_coefs_lmout(
            lm(
                exp ~ csl + csd + prl + prd + drl + rho,
                make_matrix(dfh2, 7, c( "exp", "csl", "csd", "prl", "prd", "drl", "rho"))
            ),
            out,
            c(1:6),
            dfh2$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "psl" = "Parliamentary seat share$_{t-1}$",
        "psd" = "$\\Delta$ Parliamentary seat share", 
        "csl" = "Cabinet seat share$_{t-1}$",
        "csd" = "$\\Delta$ Cabinet seat share",
        "prl" = "Relative radicalism$_{t-1}$",
        "prd" = "$\\Delta$ Relative radicalism",
        "drl" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(sel, 21, 23),
        data.frame(values = c(glance(out)$logLik, glance(out)$AIC, NA))
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccc"
    } else if (knitr::is_latex_output()){
        "lddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Selection" = 2, "Outcome" = 1), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 3)))
```

```{r}
#| label: tbl-three-ches-lr
#| tbl-cap: Three-step Heckman models accounting for absolute radicalisation based on CHES left-right radicalisation
#| results: "asis"

dfhs <- dfh |> filter(!is.na(lrgen_r_dif))
dfhc <- dfhs |> filter(cab_lag == 1)

sel <- selection(
        cab_lag ~ par_seat_share_l + lrgen_rl,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + lrgen_rl + lrgen_r_dif,
        data = dfhs
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfhs, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + lrgen_rl + lrgen_r_dif + dif_rou_l + imr,
    data = dfh2
)

modelsummary(
    list(
        "Initial Inclusion" = change_coefs(
            lm(
                exp ~ psl + prl,
                make_matrix(dfhs, 3, c("exp", "psl", "prl"))
            ),
            sel,
            c(1:3),
            type = "sel",
            cluster= dfhs$country_name
        ),
        "Continued Inclusion" = change_coefs(
            lm(
                exp ~ psl + psd + prl + prd + rho,
                make_matrix(dfhc, 6, c("exp", "psl", "psd", "prl", "prd", "rho"))
            ),
            sel,
            c(4:8),
            type = "out",
            cluster= dfhs$country_name
        ),
        "$\\Delta$ Gamson Deviation" = change_coefs_lmout(
            lm(
                exp ~ csl + csd + prl + prd + drl + rho,
                make_matrix(dfh2, 7, c( "exp", "csl", "csd", "prl", "prd", "drl", "rho"))
            ),
            out,
            c(1:6),
            dfh2$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "psl" = "Parliamentary seat share$_{t-1}$",
        "psd" = "$\\Delta$ Parliamentary seat share", 
        "csl" = "Cabinet seat share$_{t-1}$",
        "csd" = "$\\Delta$ Cabinet seat share",
        "prl" = "Absolute radicalism$_{t-1}$",
        "prd" = "$\\Delta$ Absolute radicalism",
        "drl" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(sel, 21, 23),
        data.frame(values = c(glance(out)$logLik, glance(out)$AIC, NA))
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccc"
    } else if (knitr::is_latex_output()){
        "lddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Selection" = 2, "Outcome" = 1), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 3)))
```

```{r}
#| label: tbl-three-ches-eco
#| tbl-cap: Three-step Heckman models accounting for absolute radicalisation based on CHES economic radicalisation
#| results: "asis"

dfhs <- dfh |> filter(!is.na(lrecon_r_dif))
dfhc <- dfhs |> filter(cab_lag == 1)

sel <- selection(
        cab_lag ~ par_seat_share_l + lrecon_rl,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + lrecon_rl + lrecon_r_dif,
        data = dfhs
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfhs, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + lrecon_rl + lrecon_r_dif + dif_rou_l + imr,
    data = dfh2
)

modelsummary(
    list(
        "Initial Inclusion" = change_coefs(
            lm(
                exp ~ psl + prl,
                make_matrix(dfhs, 3, c("exp", "psl", "prl"))
            ),
            sel,
            c(1:3),
            type = "sel",
            cluster= dfhs$country_name
        ),
        "Continued Inclusion" = change_coefs(
            lm(
                exp ~ psl + psd + prl + prd + rho,
                make_matrix(dfhc, 6, c("exp", "psl", "psd", "prl", "prd", "rho"))
            ),
            sel,
            c(4:8),
            type = "out",
            cluster= dfhs$country_name
        ),
        "$\\Delta$ Gamson Deviation" = change_coefs_lmout(
            lm(
                exp ~ csl + csd + prl + prd + drl + rho,
                make_matrix(dfh2, 7, c( "exp", "csl", "csd", "prl", "prd", "drl", "rho"))
            ),
            out,
            c(1:6),
            dfh2$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "psl" = "Parliamentary seat share$_{t-1}$",
        "psd" = "$\\Delta$ Parliamentary seat share", 
        "csl" = "Cabinet seat share$_{t-1}$",
        "csd" = "$\\Delta$ Cabinet seat share",
        "prl" = "Absolute radicalism$_{t-1}$",
        "prd" = "$\\Delta$ Absolute radicalism",
        "drl" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(
            sel,
            21, 23
        ),
        data.frame(
            values = c(
                glance(out)$logLik,
                glance(out)$AIC,
                NA
            )
        )
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccc"
    } else if (knitr::is_latex_output()){
        "lddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Selection" = 2, "Outcome" = 1), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 3)))
```

```{r}
#| label: tbl-three-ches-gta
#| tbl-cap: Three-step Heckman models accounting for absolute radicalisation based on CHES cultural radicalisation
#| results: "asis"

dfhs <- dfh |> filter(!is.na(galtan_r_dif))
dfhc <- dfhs |> filter(cab_lag == 1)

sel <- selection(
        cab_lag ~ par_seat_share_l + galtan_rl,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + galtan_rl + galtan_r_dif,
        data = dfhs
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfhs, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + galtan_rl + galtan_r_dif + dif_rou_l + imr,
    data = dfh2
)

modelsummary(
    list(
        "Initial Inclusion" = change_coefs(
            lm(
                exp ~ psl + prl,
                make_matrix(dfhs, 3, c("exp", "psl", "prl"))
            ),
            sel,
            c(1:3),
            type = "sel",
            cluster= dfhs$country_name
        ),
        "Continued Inclusion" = change_coefs(
            lm(
                exp ~ psl + psd + prl + prd + rho,
                make_matrix(dfhc, 6, c("exp", "psl", "psd", "prl", "prd", "rho"))
            ),
            sel,
            c(4:8),
            type = "out",
            cluster= dfhs$country_name
        ),
        "$\\Delta$ Gamson Deviation" = change_coefs_lmout(
            lm(
                exp ~ csl + csd + prl + prd + drl + rho,
                make_matrix(dfh2, 7, c( "exp", "csl", "csd", "prl", "prd", "drl", "rho"))
            ),
            out,
            c(1:6),
            dfh2$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "psl" = "Parliamentary seat share$_{t-1}$",
        "psd" = "$\\Delta$ Parliamentary seat share", 
        "csl" = "Cabinet seat share$_{t-1}$",
        "csd" = "$\\Delta$ Cabinet seat share",
        "prl" = "Absolute radicalism$_{t-1}$",
        "prd" = "$\\Delta$ Absolute radicalism",
        "drl" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(
            sel,
            21, 23
        ),
        data.frame(
            values = c(
                glance(out)$logLik,
                glance(out)$AIC,
                NA
            )
        )
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccc"
    } else if (knitr::is_latex_output()){
        "lddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Selection" = 2, "Outcome" = 1), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 3)))
```

```{r}
#| label: tbl-three-cmp-lr
#| tbl-cap: Three-step Heckman models accounting for absolute radicalisation based on CMP left-right radicalisation
#| results: "asis"

dfhs <- dfh |> filter(!is.na(rile_r_dif))
dfhc <- dfhs |> filter(cab_lag == 1)

sel <- selection(
        cab_lag ~ par_seat_share_l + rile_rl + country_name - 1,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + rile_rl + rile_r_dif,
        data = dfhs
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfhs, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + rile_rl + rile_r_dif + dif_rou_l + imr,
    data = dfh2
)

modelsummary(
    list(
        "Initial Inclusion" = change_coefs(
            lm(
                exp ~ psl + prl - 1,
                make_matrix(dfhs, 3, c("exp", "psl", "prl"))
            ),
            sel,
            c(1,2),
            type = "sel",
            cluster= dfhs$country_name
        ),
        "Continued Inclusion" = change_coefs(
            lm(
                exp ~ psl + psd + prl + prd + rho,
                make_matrix(dfhc, 6, c("exp", "psl", "psd", "prl", "prd", "rho"))
            ),
            sel,
            c(34:38),
            type = "out",
            cluster= dfhs$country_name
        ),
        "$\\Delta$ Gamson Deviation" = change_coefs_lmout(
            lm(
                exp ~ csl + csd + prl + prd + drl + rho,
                make_matrix(dfh2, 7, c( "exp", "csl", "csd", "prl", "prd", "drl", "rho"))
            ),
            out,
            c(1:6),
            dfh2$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "psl" = "Parliamentary seat share$_{t-1}$",
        "psd" = "$\\Delta$ Parliamentary seat share", 
        "csl" = "Cabinet seat share$_{t-1}$",
        "csd" = "$\\Delta$ Cabinet seat share",
        "prl" = "Absolute radicalism$_{t-1}$",
        "prd" = "$\\Delta$ Absolute radicalism",
        "drl" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(
            sel,
            21, 23
        ),
        data.frame(
            values = c(
                glance(out)$logLik,
                glance(out)$AIC,
                NA
            )
        )
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccc"
    } else if (knitr::is_latex_output()){
        "lddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Selection" = 2, "Outcome" = 1), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 3)))
```

```{r}
#| label: tbl-three-cmp-eco
#| tbl-cap: Three-step Heckman models accounting for absolute radicalisation based on CMP economic radicalisation
#| results: "asis"

dfhs <- dfh |> filter(!is.na(market_nat_r_dif))
dfhc <- dfhs |> filter(cab_lag == 1)

sel <- selection(
        cab_lag ~ par_seat_share_l + market_nat_rl + country_name - 1,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + market_nat_rl + market_nat_r_dif,
        data = dfhs
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfhs, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + market_nat_rl + market_nat_r_dif + dif_rou_l + imr,
    data = dfh2
)

modelsummary(
    list(
        "Initial Inclusion" = change_coefs(
            lm(
                exp ~ psl + prl - 1,
                make_matrix(dfhs, 3, c("exp", "psl", "prl"))
            ),
            sel,
            c(1,2),
            type = "sel",
            cluster= dfhs$country_name
        ),
        "Continued Inclusion" = change_coefs(
            lm(
                exp ~ psl + psd + prl + prd + rho,
                make_matrix(dfhc, 6, c("exp", "psl", "psd", "prl", "prd", "rho"))
            ),
            sel,
            c(34:38),
            type = "out",
            cluster= dfhs$country_name
        ),
        "$\\Delta$ Gamson Deviation" = change_coefs_lmout(
            lm(
                exp ~ csl + csd + prl + prd + drl + rho,
                make_matrix(dfh2, 7, c( "exp", "csl", "csd", "prl", "prd", "drl", "rho"))
            ),
            out,
            c(1:6),
            dfh2$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "psl" = "Parliamentary seat share$_{t-1}$",
        "psd" = "$\\Delta$ Parliamentary seat share", 
        "csl" = "Cabinet seat share$_{t-1}$",
        "csd" = "$\\Delta$ Cabinet seat share",
        "prl" = "Absolute radicalism$_{t-1}$",
        "prd" = "$\\Delta$ Absolute radicalism",
        "drl" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(
            sel,
            21, 23
        ),
        data.frame(
            values = c(
                glance(out)$logLik,
                glance(out)$AIC,
                NA
            )
        )
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccc"
    } else if (knitr::is_latex_output()){
        "lddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Selection" = 2, "Outcome" = 1), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 3)))
```

```{r}
#| label: tbl-three-cmp-gta
#| tbl-cap: Three-step Heckman models accounting for absolute radicalisation based on CMP cultural radicalisation
#| results: "asis"

dfhs <- dfh |> filter(!is.na(culture_nat_r_dif))
dfhc <- dfhs |> filter(cab_lag == 1)

sel <- selection(
        cab_lag ~ par_seat_share_l + culture_nat_rl + country_name - 1,
        as.logical(cabinet_party) ~ par_seat_share_l + par_seat_share_dif + culture_nat_rl + culture_nat_r_dif,
        data = dfhs
    )
imr <- data.frame(imr = imr_sel(sel))
dfh2 <- cbind(dfhs, imr) |>
    filter(cab_lag == 1 & cabinet_party == 1)
out <- lm(
    dif_rou_dif ~ cab_seat_share_l + cab_seat_share_dif + culture_nat_rl + culture_nat_r_dif + dif_rou_l + imr,
    data = dfh2
)

modelsummary(
    list(
        "Initial Inclusion" = change_coefs(
            lm(
                exp ~ psl + prl - 1,
                make_matrix(dfhs, 3, c("exp", "psl", "prl"))
            ),
            sel,
            c(1,2),
            type = "sel",
            cluster= dfhs$country_name
        ),
        "Continued Inclusion" = change_coefs(
            lm(
                exp ~ psl + psd + prl + prd + rho,
                make_matrix(dfhc, 6, c("exp", "psl", "psd", "prl", "prd", "rho"))
            ),
            sel,
            c(34:38),
            type = "out",
            cluster= dfhs$country_name
        ),
        "$\\Delta$ Gamson Deviation" = change_coefs_lmout(
            lm(
                exp ~ csl + csd + prl + prd + drl + rho,
                make_matrix(dfh2, 7, c( "exp", "csl", "csd", "prl", "prd", "drl", "rho"))
            ),
            out,
            c(1:6),
            dfh2$country_name
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "psl" = "Parliamentary seat share$_{t-1}$",
        "psd" = "$\\Delta$ Parliamentary seat share", 
        "csl" = "Cabinet seat share$_{t-1}$",
        "csd" = "$\\Delta$ Cabinet seat share",
        "prl" = "Absolute radicalism$_{t-1}$",
        "prd" = "$\\Delta$ Absolute radicalism",
        "drl" = "Gamson Deviation$_{t-1}$",
        "(Intercept)" = "Intercept",
        "rho" = "$\\rho$"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        heck_gof(
            sel,
            21, 23
        ),
        data.frame(
            values = c(
                glance(out)$logLik,
                glance(out)$AIC,
                NA
            )
        )
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lccc"
    } else if (knitr::is_latex_output()){
        "lddd"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Selection" = 2, "Outcome" = 1), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(16, hline_after = hl, extra_css = css) |>
row_spec(c(19:22), align = c("l", rep("c", 3)))
```

```{r}
#| label: fig-rad-gov-incl
#| fig-cap: "Marginal effect of radicalism on coalition inclusion likelihood. **Note**: 95 per cent confidence interval shown as ribbon. Marks above x-axis represent observations at according level."
#| fig-height: 3.37
#| fig-width: 6.1
#| fig-pos: t

slopes(
    lmer(cabinet_party ~ parlmean_dif_r*year + par_seat_share + (1 | country_name), df),
    datagrid(year = seq(1945,2021,1)),
    variables = "parlmean_dif_r"
) |>
mutate(type = "Relative radicalism") |>
rename(xvar = parlmean_dif_r) |>
rbind(
    slopes(
        lmer(cabinet_party ~ left_right_r*year + par_seat_share + (1 | country_name), df),
        datagrid(year = seq(1945,2021,1)),
        variables = "left_right_r"
    ) |>
    mutate(type = "Absolute radicalism") |>
    rename(xvar = left_right_r)
) |>
ggplot(aes(year, estimate)) +
    facet_wrap(~type) +
    geom_hline(yintercept = 0, color = "lightcoral") +
    geom_line() +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3) +
    geom_rug(
        data = df |> select(year),
        aes(x = year, 0),
        position = position_jitter(),
        alpha = .2,
        sides = "b"
    ) +
    scale_y_continuous(
        name = "Marginal Effect",
        limits = c(-.2, 0)
    ) +
    scale_x_continuous(
        name = "",
        breaks = seq(1945, 2021, 15),
        expand = c(0,0)
    ) +
    theme(
        panel.spacing = unit(2, "lines")
    )
```

```{r}
#| label: value-data

dfv <- rbind(
    df |>
        filter(
            !is.na(highs),
            !is.na(meds),
            !is.na(lows),
            cabinet_party == 1
        ) |>
        uncount(highs) |>
        mutate(value = "high") |>
        select(-c(meds, lows)),
    df |>
        filter(
            !is.na(highs),
            !is.na(meds),
            !is.na(lows),
            cabinet_party == 1
        ) |>
        uncount(meds) |>
        mutate(value = "med") |>
        select(-c(highs, lows)),
    df |>
        filter(
            !is.na(highs),
            !is.na(meds),
            !is.na(lows),
            cabinet_party == 1
        ) |>
        uncount(lows) |>
        mutate(value = "low") |>
        select(-c(highs, meds))
) |>
mutate(value = factor(value, levels = c("low", "med", "high")))
```

```{r}
#| label: tbl-value
#| tbl-cap: Ordinal logistic model for portfolio value
#| results: "asis"

mv1 <- clm(value ~ cab_seat_share + parlmean_dif_r + dif_rou, data = dfv)
mv2 <- clm(value ~ cab_seat_share + left_right_r + dif_rou, data = dfv)

modelsummary(
    list(
        "Model 1" = coeftest(mv1, vcovCL, cluster = dfv$country_name),
        "Model 2" = coeftest(mv2, vcovCL, cluster = dfv$country_name)
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "cab_seat_share" = "Cabinet seat share",
        "parlmean_dif_r" = "Relative radicalism",
        "left_right_r" = "Absolute radicalism",
        "dif_rou" = "Gamson Deviation",
        "low|med" = "First threshold",
        "med|high" = "Second threshold"
    ),
    gof_map = list(
        list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f"),
        list("raw" = "aic", "clean" = "AIC", "fmt" = "%.1f"),
        list("raw" = "bic", "clean" = "BIC", "fmt" = "%.1f")
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$; Country robust standard errors in parantheses."),
    align = if (knitr::is_html_output()){
        "lcc"
    } else if (knitr::is_latex_output()){
        "ldd"
    },
    output = "latex",
    escape = FALSE
) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(8, hline_after = hl, extra_css = css) |>
row_spec(c(13:15), align = c("l", rep("c", 2)))
```

```{r}
#| label: fig-value
#| fig-cap: "Predicted probability of portfolio value. **Note:** Predicted probabilities based on ordinal logistic models. 95 per cent confidence interval shown as ribbon."
#| fig-height: 4
#| fig-width: 6.1

rbind(
    predictions(
        mv1,
        newdata = datagrid(parlmean_dif_r = seq(0,5,.1)),
        vcov = vcovCL(mv1, cluster = dfv$country_name)
    ) |>
    rename(xvar = parlmean_dif_r) |>
    mutate(
        var = "Radicalism",
        model = "Model 1 (Radicalism = Relative)"
    ) |>
    select(-c(dif_rou, cab_seat_share)),
    predictions(
        mv1, 
        newdata = datagrid(cab_seat_share = seq(0,1,.01)),
        vcov = vcovCL(mv1, cluster = dfv$country_name)
    ) |>
    rename(xvar = cab_seat_share) |>
    mutate(
        var = "Cabinet Seat Share",
        model = "Model 1 (Radicalism = Relative)"
    ) |>
    select(-c(parlmean_dif_r, dif_rou)),
    predictions(
        mv1,
        newdata = datagrid(dif_rou = seq(-50,50,1)),
        vcov = vcovCL(mv1, cluster = dfv$country_name)
    ) |>
    rename(xvar = dif_rou) |>
    mutate(
        var = "Gamson Deviation",
        model = "Model 1 (Radicalism = Relative)"
    ) |>
    select(-c(parlmean_dif_r, cab_seat_share)),
    predictions(
        mv2,
        newdata = datagrid(left_right_r = seq(0,5,.1)),
        vcov = vcovCL(mv2, cluster = dfv$country_name)
    ) |>
    rename(xvar = left_right_r) |>
    mutate(
        var = "Radicalism",
        model = "Model 2 (Radicalism = Absolute)"
    ) |>
    select(-c(dif_rou, cab_seat_share)),
    predictions(
        mv2,
        newdata = datagrid(cab_seat_share = seq(0,1,.01)),
        vcov = vcovCL(mv2, cluster = dfv$country_name)
    ) |>
    rename(xvar = cab_seat_share) |>
    mutate(
        var = "Cabinet Seat Share",
        model = "Model 2 (Radicalism = Absolute)"
    ) |>
    select(-c(left_right_r, dif_rou)),
    predictions(
        mv2,
        newdata = datagrid(dif_rou = seq(-50,50,1)),
        vcov = vcovCL(mv2, cluster = dfv$country_name)
    ) |>
    rename(xvar = dif_rou) |>
    mutate(
        var = "Gamson Deviation",
        model = "Model 2 (Radicalism = Absolute)"
    ) |>
    select(-c(left_right_r, cab_seat_share))
) |>
mutate(model = str_wrap(model, width = 20)) |>
ggplot(aes(xvar, estimate, colour = group, fill = group)) +
    facet_grid(model ~ var, scales = "free_x") +
    geom_line() +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), colour = NA, alpha = .3) +
    scale_y_continuous(
        name = "Predicted Probability",
        limits = c(0,.75),
        breaks = seq(0,.75,.25),
        expand = c(0,0),
        labels = scales::label_percent()
    ) +
    scale_x_continuous(
        name = "",
        expand = c(0,0)
    ) +
    scale_fill_discrete(name = "Portfolio Value") +
    scale_colour_discrete(name = "Portfolio Value") +
    theme(
        panel.spacing = unit(2, "lines"),
        legend.position = "bottom",
    )
```

```{r}
#| label: fig-zero-infl
#| fig-cap: "Distribution rounded dependent variable. **Note**: Binwidth is set to 1."
#| fig-height: 3.37
#| fig-width: 6.1

dfc |>
    ggplot(aes(dif_rou)) +
        geom_histogram(binwidth = 1) +
        scale_x_continuous(
            name = "Gamson Deviation",
            expand = c(0,0)
        ) +
        scale_y_continuous(
            name = "Observations",
            limits = c(0,650),
            expand = c(0, 0)
        )
```

```{r}
#| label: tbl-zero-infl
#| tbl-cap: Multi-level zero-inflated models
#| results: "asis"

dfzi <- dfc |>
    group_by(cabinet_id) |>
    mutate(
        seat_dif = max(cab_seat_share) - min(cab_seat_share),
        rel_dif = max(parlmean_dif_r) - min(parlmean_dif_r),
        abs_dif = max(left_right_r) - min(left_right_r),
        cabinet_id = as.factor(cabinet_id)
    )

mzr <- glmmTMB(
    dif_rou ~ parlmean_dif_r + cab_seat_share + (0 + cab_seat_share | country_name),
    ziformula = ~ poly(parlmean_dif_r, 2, raw = TRUE) + poly(cab_seat_share, 2, raw = TRUE) + rel_dif + seat_dif + (1 | cabinet_id) + (1 + seat_dif | country_name),
    data = dfzi
)
mza <- glmmTMB(
    dif_rou ~ left_right_r + cab_seat_share + (0 + cab_seat_share | country_name),
    ziformula = ~ poly(left_right_r, 2, raw = TRUE) + poly(cab_seat_share, 2, raw = TRUE) + abs_dif + seat_dif + (1 | cabinet_id) + (1 + seat_dif | country_name),
    data = dfzi
)

modelsummary(
    list(
        "Zero-inflated portion" = change_zero(
            lm(
                zer ~ pr1 + pr2 + css1 + css2 + rd + sd,
                make_matrix(
                    dfc,
                    7,
                    c("zer", "pr1", "pr2", "css1", "css2", "rd", "sd")
                )
            ),
            mzr,
            type = "zero"
        ),
        "Conditional portion" = change_zero(
            lm(
                dif ~ pr1 + css1,
                make_matrix(
                    dfc,
                    3,
                    c("dif", "pr1", "css1")
                )
            ),
            mzr,
            type = "conditional"
        ),
        "Zero-inflated portion" = change_zero(
            lm(
                zer ~ lr1 + lr2 + css1 + css2 + ad + sd,
                make_matrix(
                    dfc,
                    7,
                    c("zer", "lr1", "lr2", "css1", "css2", "ad", "sd")
                )
            ),
            mza,
            type = "zero"
        ),
        "Conditional portion" = change_zero(
            lm(
                dif ~ lr1 + css1,
                make_matrix(
                    dfc,
                    3,
                    c("dif", "lr1", "css1")
                )
            ),
            mza,
            type = "conditional"
        )
    ),
    estimate = "{estimate}{stars}",
    stars = c("$^{*}$" = .05, "$^{**}$" = .01, "$^{***}$" = .001),
    coef_map = c(
        "pr1" = "Relative radicalism",
        "pr2" = "Relative radicalism²",
        "lr1" = "Absolute radicalism",
        "lr2" = "Absolute radicalism²",
        "css1" = "Cabinet seat share",
        "css2" = "Cabinet seat share²",
        "rd" = "Max dif. of relative radicalism",
        "ad" = "Max dif. of absolute radicalism",
        "sd" = "Max dif. of cabinet seat share",
        "(Intercept)" = "Intercept"
    ),
    gof_map = list(list("raw" = "nobs", "clean" = "Observations", "fmt" = "%.0f")),
    add_rows = cbind(
        zero_gof(mzr),
        zero_gof(mza) |> select(-rowname)
    ),
    notes = list("$^{*} p<.05$; $^{**} p<.01$; $^{***} p<.001$"),
    align = if(knitr::is_latex_output()){
        "ldddd"
    } else if(knitr::is_html_output()){
        "lcccc"
    },
    output = "latex",
    escape = FALSE
) |>
add_header_above(c("", "Model 1" = 2, "Model 2" = 2), bold = TRUE) |>
kable_styling(latex_options = ltop, full_width = FALSE) |>
row_spec(0, bold = TRUE, align = "c") |>
row_spec(c(21:26), align = c("l", rep("c", 4)))
```

```{r}
#| label: fig-zero-infl-pred
#| fig-cap: "Predictions from zero-inflation models. **Note**: Predictions from zero-inflation models. Left-hand side panels depict estimated probabilities of a Gamsonian allocation over different IV-levels. Right-hand side panels depict estimated levels of deviations from a proportional ministry allocation over these IVs. 95 per cent confidence interval depicted as ribbon. Marks on bottom and top of each panel represent observations at according level. These marks are filtered for non-zero values for the conditional part."
#| fig-height: 8
#| fig-width: 6.1

predictions(
    mzr,
    newdata = datagrid(rel_dif = seq(0,5,.1)),
    type = "zprob"
) |>
rename(xval = rel_dif) |>
mutate(
    rad = "Model 1: Relative radicalism",
    xvar = "Radicalism difference",
    mod = "Zero-inflated"
) |>
select(estimate, conf.low, conf.high, xval, rad, xvar, mod) |>
rbind(
    predictions(
        mza,
        newdata = datagrid(abs_dif = seq(0,5,.1)),
        type = "zprob"
    ) |>
    rename(xval = abs_dif) |>
    mutate(
        rad = "Model 2: Absolute radicalism",
        xvar = "Radicalism difference",
        mod = "Zero-inflated"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mzr,
        newdata = datagrid(seat_dif = seq(0,1,.01)),
        type = "zprob"
    ) |>
    rename(xval = seat_dif) |>
    mutate(
        rad = "Model 1: Relative radicalism",
        xvar = "Seat difference",
        mod = "Zero-inflated"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mza,
        newdata = datagrid(seat_dif = seq(0,1,.01)),
        type = "zprob"
    ) |>
    rename(xval = seat_dif) |>
    mutate(
        rad = "Model 2: Absolute radicalism",
        xvar = "Seat difference",
        mod = "Zero-inflated"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mzr,
        newdata = datagrid(cab_seat_share = seq(0,1,.01)),
        type = "zprob"
    ) |>
    rename(xval = cab_seat_share) |>
    mutate(
        rad = "Model 1: Relative radicalism",
        xvar = "Seat share",
        mod = "Zero-inflated"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mza,
        newdata = datagrid(cab_seat_share = seq(0,1,.01)),
        type = "zprob"
    ) |>
    rename(xval = cab_seat_share) |>
    mutate(
        rad = "Model 2: Absolute radicalism",
        xvar = "Seat share",
        mod = "Zero-inflated"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mzr,
        newdata = datagrid(cab_seat_share = seq(0,1,.01)),
        type = "conditional"
    ) |>
    rename(xval = cab_seat_share) |>
    mutate(
        rad = "Model 1: Relative radicalism",
        xvar = "Seat share",
        mod = "Conditional"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mza,
        newdata = datagrid(cab_seat_share = seq(0,1,.01)),
        type = "conditional"
    ) |>
    rename(xval = cab_seat_share) |>
    mutate(
        rad = "Model 2: Absolute radicalism",
        xvar = "Seat share",
        mod = "Conditional"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
)  |>
rbind(
    predictions(
        mzr,
        newdata = datagrid(parlmean_dif_r = seq(0,5,.1)),
        type = "zprob"
    ) |>
    rename(xval = parlmean_dif_r) |>
    mutate(
        rad = "Model 1: Relative radicalism",
        xvar = "Radicalism",
        mod = "Zero-inflated"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mza,
        newdata = datagrid(left_right_r = seq(0,5,.1)),
        type = "zprob"
    ) |>
    rename(xval = left_right_r) |>
    mutate(
        rad = "Model 2: Absolute radicalism",
        xvar = "Radicalism",
        mod = "Zero-inflated"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mzr,
        newdata = datagrid(parlmean_dif_r = seq(0,5,.1)),
        type = "conditional"
    ) |>
    rename(xval = parlmean_dif_r) |>
    mutate(
        rad = "Model 1: Relative radicalism",
        xvar = "Radicalism",
        mod = "Conditional"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
rbind(
    predictions(
        mza,
        newdata = datagrid(left_right_r = seq(0,5,.1)),
        type = "conditional"
    ) |>
    rename(xval = left_right_r) |>
    mutate(
        rad = "Model 2: Absolute radicalism",
        xvar = "Radicalism",
        mod = "Conditional"
    ) |>
    select(estimate, conf.low, conf.high, xval, rad, xvar, mod)
) |>
mutate(
    mod = factor(
        ifelse(
            mod == "Zero-inflated",
            "Zero-inflated\nPredicted probability of\nproportional allocation",
            "Conditional\nPredicted Gamson deviation"
        ),            
        levels = c(
            "Zero-inflated\nPredicted probability of\nproportional allocation",
            "Conditional\nPredicted Gamson deviation"
        )
    ),
    xvar = factor(
        xvar,
        levels = c(
            "Radicalism difference",
            "Seat difference",
            "Radicalism",
            "Seat share"
        )
    )
) |>
ggplot(aes(xval, estimate, colour = rad, fill = rad)) +
    facet_grid2(
        rows = vars(xvar), 
        cols = vars(mod),
        scales = "free",
        independent = "all"
    ) +
    geom_line() +
    geom_ribbon(
        aes(ymin = conf.low, ymax = conf.high), 
        alpha = .3,
        colour = NA
    ) +
    geom_rug(
        data = dfzi |>
            ungroup() |>
            select(rel_dif) |>
            rename(xval = rel_dif) |>
            mutate(
                xvar = "Radicalism difference",
                mod = "Zero-inflated"
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    select(seat_dif) |>
                    rename(xval = seat_dif) |>
                    mutate(
                        xvar = "Seat difference",
                        mod = "Zero-inflated"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    select(cab_seat_share) |>
                    rename(xval = cab_seat_share) |>
                    mutate(
                        xvar = "Seat share",
                        mod = "Zero-inflated"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    filter(dif_rou != 0) |>
                    select(cab_seat_share) |>
                    rename(xval = cab_seat_share) |>
                    mutate(
                        xvar = "Seat share",
                        mod = "Conditional"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    select(parlmean_dif_r) |>
                    rename(xval = parlmean_dif_r) |>
                    mutate(
                        xvar = "Radicalism",
                        mod = "Zero-inflated"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    filter(dif_rou != 0) |>
                    select(parlmean_dif_r) |>
                    rename(xval = parlmean_dif_r) |>
                    mutate(
                        xvar = "Radicalism",
                        mod = "Conditional"
                    )
            ) |>
            mutate(
                rad = "Model 1: Relative radicalism",
                mod = factor(
                    ifelse(
                        mod == "Zero-inflated",
                        "Zero-inflated\nPredicted probability of\nproportional allocation",
                        "Conditional\nPredicted Gamson deviation"
                    ),            
                    levels = c(
                        "Zero-inflated\nPredicted probability of\nproportional allocation",
                        "Conditional\nPredicted Gamson deviation"
                    )
                ),
                xvar = factor(
                    xvar,
                    levels = c(
                        "Radicalism difference",
                        "Seat difference",
                        "Radicalism",
                        "Seat share"
                    )
                )
            ), 
        y = 0, 
        sides = "b", 
        alpha = .3
    ) +
    geom_rug(
        data = dfzi |>
            ungroup() |>
            select(abs_dif) |>
            rename(xval = abs_dif) |>
            mutate(
                xvar = "Radicalism difference",
                mod = "Zero-inflated"
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    select(seat_dif) |>
                    rename(xval = seat_dif) |>
                    mutate(
                        xvar = "Seat difference",
                        mod = "Zero-inflated"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    select(cab_seat_share) |>
                    rename(xval = cab_seat_share) |>
                    mutate(
                        xvar = "Seat share",
                        mod = "Zero-inflated"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    filter(dif_rou != 0) |>
                    select(cab_seat_share) |>
                    rename(xval = cab_seat_share) |>
                    mutate(
                        xvar = "Seat share",
                        mod = "Conditional"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    select(left_right_r) |>
                    rename(xval = left_right_r) |>
                    mutate(
                        xvar = "Radicalism",
                        mod = "Zero-inflated"
                    )
            ) |>
            rbind(
                dfzi |>
                    ungroup() |>
                    filter(dif_rou != 0) |>
                    select(left_right_r) |>
                    rename(xval = left_right_r) |>
                    mutate(
                        xvar = "Radicalism",
                        mod = "Conditional"
                    )
            ) |>
            mutate(
                rad = "Model 2: Absolute radicalism",
                mod = factor(
                    ifelse(
                        mod == "Zero-inflated",
                        "Zero-inflated\nPredicted probability of\nproportional allocation",
                        "Conditional\nPredicted Gamson deviation"
                    ),            
                    levels = c(
                        "Zero-inflated\nPredicted probability of\nproportional allocation",
                        "Conditional\nPredicted Gamson deviation"
                    )
                ),
                xvar = factor(
                    xvar,
                    levels = c(
                        "Radicalism difference",
                        "Seat difference",
                        "Radicalism",
                        "Seat share"
                    )
                )
            ), 
        y = 0, 
        sides = "t", 
        alpha = .3
    ) +
    facetted_pos_scales(
        x = list(
            scale_x_continuous(
                name = "",
                limits = c(0,5),
                expand = c(0,0)
            ),
            NULL,
            scale_x_continuous(
                limits = c(0,1),
                expand = c(0,0)
            ),
            NULL,
            scale_x_continuous(
                limits = c(0,5),
                expand = c(0,0)
            ),
            scale_x_continuous(
                limits = c(0,5),
                expand = c(0,0)
            ),
            scale_x_continuous(
                limits = c(0,1),
                expand = c(0,0),
                labels = scales::label_percent()
            ),
            scale_x_continuous(
                limits = c(0,1),
                expand = c(0,0),
                labels = scales::label_percent()
            )
        ),
        y = list(
            scale_y_continuous(
                name = "",
                limits = c(0,1),
                expand = c(0,0),
                labels = scales::label_percent()
            ),
            NULL,
            scale_y_continuous(
                limits = c(0,1),
                expand = c(0,0),
                labels = scales::label_percent()
            ),
            NULL,
            scale_y_continuous(
                limits = c(0,1),
                expand = c(0,0),
                labels = scales::label_percent()
            ),
            scale_y_continuous(
                limits = c(-20,10),
                expand = c(0,0)
            ),
            scale_y_continuous(
                limits = c(0,1),
                expand = c(0,0),
                labels = scales::label_percent()
            ),
            scale_y_continuous(
                limits = c(-20,10),
                expand = c(0,0)
            )
        )
    ) +
    theme(
        legend.position = "bottom",
        legend.title = element_blank(),
        panel.spacing = unit(1.25, "lines")
    )
```